cgma
Public Member Functions | Private Attributes
CubitMatrix Class Reference

#include <CubitMatrix.hpp>

Inheritance diagram for CubitMatrix:
CubitTransformMatrix

List of all members.

Public Member Functions

 CubitMatrix ()
 CubitMatrix (const int num_rows, const int num_cols)
 CubitMatrix (const CubitVector &vec1, const CubitVector &vec2, const CubitVector &vec3)
 CubitMatrix (const CubitVector &vec1, const CubitVector &vec2, const CubitVector &vec3, const CubitVector &vec4)
 CubitMatrix (const CubitVector &vec1, const CubitVector &vec2)
 CubitMatrix (const int n)
 CubitMatrix (std::vector< int > &is, std::vector< int > &js, std::vector< double > &es, int n, int m)
 CubitMatrix (const CubitMatrix &matrix)
virtual ~CubitMatrix ()
void reset_size (const int n, const int m, double default_value=0.0)
int num_rows () const
int num_cols () const
void print_matrix () const
void print_matrix (char *filename) const
void set_to_identity ()
bool is_identity (double tol=1e-8) const
bool is_equal (const CubitMatrix &other, double tol=1e-8) const
void plus_identity ()
void set (const int n, const int m, const double val)
void add (const int n, const int m, const double val)
double get (int n, int m) const
CubitMatrix operator= (const CubitMatrix &matrix)
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)
CubitMatrix inverse ()
CubitMatrix symm () const
CubitBoolean positive_definite () const
double determinant () const
double inf_norm () const
double frobenius_norm_squared () const
double frobenius_norm_squared_symm () const
double frobenius_norm_squared_adj () const
double frobenius_norm_squared_inv () const
double condition () const
int rank (void) const
double cofactor (const int row, const int col) const
CubitMatrix adjoint () const
CubitMatrix transpose () const
CubitMatrix sub_matrix (const int row, const int col) const
void sub_matrix (const std::vector< bool > &rows_to_include, const std::vector< bool > &cols_to_include, CubitMatrix &submatrix)
int gauss_elim (CubitVector &b)
int factor (CubitVector &pivot)
void solve (CubitVector &b, const CubitVector &pivot)
CubitStatus solveNxN (CubitMatrix &rhs, CubitMatrix &coef)
CubitStatus solveNxN (const std::vector< double > &rhs, std::vector< double > &coef)
CubitStatus ludcmp (double *indx, double &d)
CubitStatus lubksb (double *indx, double *b)
bool vector_outer_product (const CubitVector &vec1, const CubitVector &vec2)

Private Attributes

double ** matrixPtr
double * matrixMem
int numRows
int numCols

Detailed Description

Definition at line 20 of file CubitMatrix.hpp.

Constructor & Destructor Documentation

CubitMatrix::CubitMatrix ( )

Definition at line 14 of file CubitMatrix.cpp.

{
  matrixMem = NULL;
  matrixPtr = NULL;
  reset_size( 3, 3 );
}
CubitMatrix::CubitMatrix ( const int  num_rows,
const int  num_cols 
)

Definition at line 21 of file CubitMatrix.cpp.

{
  matrixMem = NULL;
  matrixPtr = NULL;
  reset_size( n, m );
}
CubitMatrix::CubitMatrix ( const CubitVector &  vec1,
const CubitVector &  vec2,
const CubitVector &  vec3 
)

Definition at line 39 of file CubitMatrix.cpp.

{
  matrixMem = NULL;
  matrixPtr = NULL;
  reset_size( 3, 3 );
  
    // Initialize the matrix columns to the three vectors
  matrixPtr[0][0] = vec1.x();
  matrixPtr[1][0] = vec1.y();
  matrixPtr[2][0] = vec1.z();
  matrixPtr[0][1] = vec2.x();
  matrixPtr[1][1] = vec2.y();
  matrixPtr[2][1] = vec2.z();
  matrixPtr[0][2] = vec3.x();
  matrixPtr[1][2] = vec3.y();
  matrixPtr[2][2] = vec3.z();
}
CubitMatrix::CubitMatrix ( const CubitVector &  vec1,
const CubitVector &  vec2,
const CubitVector &  vec3,
const CubitVector &  vec4 
)

Definition at line 59 of file CubitMatrix.cpp.

{
  matrixMem = NULL;
  matrixPtr = NULL;
  reset_size( 3, 4 );
  
    // Initialize the matrix columns to the four vectors
  matrixPtr[0][0] = vec1.x();
  matrixPtr[1][0] = vec1.y();
  matrixPtr[2][0] = vec1.z();
  matrixPtr[0][1] = vec2.x();
  matrixPtr[1][1] = vec2.y();
  matrixPtr[2][1] = vec2.z();
  matrixPtr[0][2] = vec3.x();
  matrixPtr[1][2] = vec3.y();
  matrixPtr[2][2] = vec3.z();
  matrixPtr[0][3] = vec4.x();
  matrixPtr[1][3] = vec4.y();
  matrixPtr[2][3] = vec4.z();
}
CubitMatrix::CubitMatrix ( const CubitVector &  vec1,
const CubitVector &  vec2 
)

Definition at line 83 of file CubitMatrix.cpp.

{             
  matrixMem = NULL;
  matrixPtr = NULL;
  reset_size( 3, 3 );
 
  this->vector_outer_product(vec1, vec2); 
}
CubitMatrix::CubitMatrix ( const int  n)

Definition at line 29 of file CubitMatrix.cpp.

{
  matrixMem = NULL;
  matrixPtr = NULL;
  reset_size( n, n );
  
    // Initialize matrix to identity.
  set_to_identity();
}
CubitMatrix::CubitMatrix ( std::vector< int > &  is,
std::vector< int > &  js,
std::vector< double > &  es,
int  n,
int  m 
)

Definition at line 94 of file CubitMatrix.cpp.

{
  matrixMem = NULL;
  matrixPtr = NULL;
  reset_size( n, m );
  int length = is.size();
  for ( int k = 0; k < length; k++ )
  {
    int i = is[k];
    int j = js[k];
    matrixPtr[i][j] += es[k];
  }
}
CubitMatrix::CubitMatrix ( const CubitMatrix &  matrix)

Definition at line 114 of file CubitMatrix.cpp.

{
  matrixMem = NULL;
  matrixPtr = NULL;
  reset_size( matrix.num_rows(), matrix.num_cols() );

  int ii;
  for( ii = 0; ii < numRows; ii++ )
  {
    for( int jj = 0; jj < numCols; jj++ )
       matrixPtr[ii][jj] = matrix.get( ii, jj );
  }
}
CubitMatrix::~CubitMatrix ( ) [virtual]

Definition at line 130 of file CubitMatrix.cpp.

{
  delete [] matrixPtr;
  delete [] matrixMem;
}

Member Function Documentation

void CubitMatrix::add ( const int  n,
const int  m,
const double  val 
) [inline]

Definition at line 91 of file CubitMatrix.hpp.

    {
      matrixPtr[n][m] += val;
    }
CubitMatrix CubitMatrix::adjoint ( ) const

Definition at line 566 of file CubitMatrix.cpp.

{
  CubitMatrix matrix( numRows, numRows );
  
  for( int ii = 0; ii < numRows; ii++ )
  {
    for( int jj = 0; jj < numRows; jj++ )
    {
      matrix.set( ii, jj, cofactor( ii, jj ) );
    }
  }
  return matrix.transpose();
}
double CubitMatrix::cofactor ( const int  row,
const int  col 
) const

Definition at line 553 of file CubitMatrix.cpp.

{
  double c = 0.0;
  CubitMatrix matrix_sub( numRows - 1, numCols -1 );
  
  matrix_sub = sub_matrix( row, col );
  
  c = matrix_sub.determinant();
  c = (row+col)%2 ? -1*c : c;
  
  return c;
}
double CubitMatrix::condition ( ) const

Definition at line 756 of file CubitMatrix.cpp.

{ 
    // condition number of A using frobenius norm 
  
  double norm = ( this->frobenius_norm_squared() ) * (this->frobenius_norm_squared_inv() );

  return sqrt( norm );

}
double CubitMatrix::determinant ( ) const

Definition at line 527 of file CubitMatrix.cpp.

{
  double det = 0.0;
  
  if( numRows == 1 )
     det = matrixPtr[0][0];
  else if( numRows == 2 )
     det = matrixPtr[0][0] * matrixPtr[1][1] - matrixPtr[0][1]
        * matrixPtr[1][0];
  else if (numRows == 3)
     det = matrixPtr[0][0] * matrixPtr[1][1] * matrixPtr[2][2] +
         matrixPtr[0][1] * matrixPtr[1][2] * matrixPtr[2][0] +
         matrixPtr[0][2] * matrixPtr[1][0] * matrixPtr[2][1] -
         matrixPtr[2][0] * matrixPtr[1][1] * matrixPtr[0][2] -
         matrixPtr[2][1] * matrixPtr[1][2] * matrixPtr[0][0] -
         matrixPtr[2][2] * matrixPtr[1][0] * matrixPtr[0][1];  
  else
  {
    for( int jj = 0; jj < numRows; jj++ )
    {
      det += ( matrixPtr[0][jj] * cofactor( 0, jj ) );
    }
  }
  return det;
}
int CubitMatrix::factor ( CubitVector &  pivot)

Definition at line 835 of file CubitMatrix.cpp.

{
  double pvt[3];

  const int n=3;
  double s[3], tmp;

  int i,j;
  for ( i=0; i<n; i++ ) 
  {
    s[i] = 0.0;
    for ( j=0; j<n; j++ ) 
    {
      tmp = fabs( matrixPtr[i][j] );
      if ( tmp > s[i] ) 
      {
        s[i] = tmp;
      }
    }
       
    if ( s[i] == 0.0 )
    {
      return(1);
    }
  }

  for ( int k=0; k<n-1; k++ ) 
  {
    double ck = 0.0;
    int i0 = -1;
    for ( i=k; i<n; i++ ) 
    {
      tmp = fabs( matrixPtr[i][k] / s[i] );
      if ( tmp > ck ) 
      {
        ck = tmp;
        i0 = i;
      }
    }

    pvt[k] = i0;
    if ( ck == 0.0 ) { return(1); }

    if ( i0 != k ) 
    {
      for ( j=k; j<n; j++ ) 
      {
        double swap = matrixPtr[i0][j];
        matrixPtr[i0][j] = matrixPtr[k][j];
        matrixPtr[k][j] = swap;
      }
    }

    for ( i=k+1; i<n; i++ ) 
    {
      double r = matrixPtr[i][k] / matrixPtr[k][k];
      matrixPtr[i][k] = r;
      for ( j=k+1; j<n; j++ ) 
      {
        matrixPtr[i][j] -= r * matrixPtr[k][j];
      }
    }
  }

  pivot.set( pvt[0], pvt[1], pvt[2] );
  return(0);
}
double CubitMatrix::frobenius_norm_squared ( ) const

Definition at line 674 of file CubitMatrix.cpp.

{ 
    // frobenius norm-squared  = trace( M^T M )
  
  double matrix_norm=0;
  for ( int ii = 0; ii < numRows; ii++ ) {

    for( int jj = 0; jj < numCols; jj++ )
    {
      matrix_norm += matrixPtr[ii][jj] * matrixPtr[ii][jj];
    }

  }

  return matrix_norm;

}
double CubitMatrix::frobenius_norm_squared_adj ( ) const

Definition at line 715 of file CubitMatrix.cpp.

{ 
    // square of frobenius norm of adjoint
  
  double norm=0;

  if ( numRows == 1 ) { norm=1; }

  if ( numRows == 2 ) {
    norm = this->frobenius_norm_squared();
  }

  if ( numRows == 3 ) {
    norm = 0.5 * ( pow( this->frobenius_norm_squared(), 2 ) - this->frobenius_norm_squared_symm() );
  }

  if ( numRows > 3 ) {
    CubitMatrix adj = this->adjoint();
    norm = adj.frobenius_norm_squared();
  }

  return norm;

}
double CubitMatrix::frobenius_norm_squared_inv ( ) const

Definition at line 740 of file CubitMatrix.cpp.

{ 
    // square of frobenius norm of A-inverse
  
  double det = this->determinant();
  
  if(det==0)
    throw std::invalid_argument ("Determinant cannot be 0");
  //assert( det != 0 );

  double norm=this->frobenius_norm_squared_adj()/pow(det,2);

  return norm;

}
double CubitMatrix::frobenius_norm_squared_symm ( ) const

Definition at line 692 of file CubitMatrix.cpp.

{ 
    // frobenius norm-squared 2 = trace[( M^T M )( M^T M )]
  
  double matrix_norm=0;
  for ( int ii = 0; ii < numRows; ii++ )
  {
    for( int jj = 0; jj < numCols; jj++ )
    {
      double b=0;
      for ( int kk = 0; kk < numRows; kk++ )
      {
        b += matrixPtr[kk][ii] * matrixPtr[kk][jj];
      }
      matrix_norm += b*b;
    }

  }

  return matrix_norm;

}
int CubitMatrix::gauss_elim ( CubitVector &  b)

Definition at line 827 of file CubitMatrix.cpp.

{
    CubitVector pivot;
    int ierr = factor( pivot );
    if ( ierr == 0 ) {  solve( b, pivot ); }
    return ierr;
}
double CubitMatrix::get ( int  n,
int  m 
) const [inline]

Definition at line 97 of file CubitMatrix.hpp.

    {
      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];
    }
double CubitMatrix::inf_norm ( ) const

Definition at line 657 of file CubitMatrix.cpp.

{ 
    // infinity norm  = max_i sum_j  | A_ij |
  double matrix_norm = 0., row_norm, v;
  for ( int ii = 0; ii < numRows; ii++ ) {
    row_norm = 0.;
    for( int jj = 0; jj < numCols; jj++ )
    {
      v = fabs( get( ii, jj ) );
      row_norm += v;
    }
    if ( row_norm > matrix_norm )
       matrix_norm = row_norm;
  }
  return matrix_norm;
}
CubitMatrix CubitMatrix::inverse ( )

Reimplemented in CubitTransformMatrix.

Definition at line 429 of file CubitMatrix.cpp.

{
  // can't invert a non-square matrix
  if(numRows!=numCols)
    throw std::invalid_argument ("Cannot invert a non-Square matrix");
  //assert(numRows == numCols);

  CubitMatrix matrix_inverse( numRows, numCols );
  
  if (numRows <4)
  {
    double   det;
    det = determinant();
    if(fabs(det) <= CUBIT_DBL_MIN)
      throw std::invalid_argument ("Determinants Absolute value must be greater that CUBIT_DBL_MIN");
    //assert( fabs(det) > CUBIT_DBL_MIN );
    double det_inv = 1./det;

    if ( numRows == 1 ) {
      det = determinant();
      if(fabs(det) <= CUBIT_DBL_MIN)
      throw std::invalid_argument ("Determinants Absolute value must be greater that CUBIT_DBL_MIN");
      //assert( fabs(det) > CUBIT_DBL_MIN );
    
      matrix_inverse.set(0,0, matrixPtr[0][0]);
    }

    if ( numRows == 2 ) {
      matrix_inverse.set(0,0, matrixPtr[1][1]);
      matrix_inverse.set(1,0,-matrixPtr[1][0]);
      matrix_inverse.set(0,1,-matrixPtr[0][1]);
      matrix_inverse.set(1,1, matrixPtr[0][0]);
    }

    if ( numRows == 3 ) {
      matrix_inverse.set(0,0, matrixPtr[1][1] * matrixPtr[2][2] - matrixPtr[1][2] * matrixPtr[2][1] );
      matrix_inverse.set(1,0, matrixPtr[2][0] * matrixPtr[1][2] - matrixPtr[1][0] * matrixPtr[2][2] );
      matrix_inverse.set(2,0, matrixPtr[1][0] * matrixPtr[2][1] - matrixPtr[1][1] * matrixPtr[2][0] );
      matrix_inverse.set(0,1, matrixPtr[2][1] * matrixPtr[0][2] - matrixPtr[0][1] * matrixPtr[2][2] );
      matrix_inverse.set(1,1, matrixPtr[0][0] * matrixPtr[2][2] - matrixPtr[0][2] * matrixPtr[2][0] );
      matrix_inverse.set(2,1, matrixPtr[0][1] * matrixPtr[2][0] - matrixPtr[0][0] * matrixPtr[2][1] );
      matrix_inverse.set(0,2, matrixPtr[0][1] * matrixPtr[1][2] - matrixPtr[0][2] * matrixPtr[1][1] );
      matrix_inverse.set(1,2, matrixPtr[1][0] * matrixPtr[0][2] - matrixPtr[0][0] * matrixPtr[1][2] );
      matrix_inverse.set(2,2, matrixPtr[0][0] * matrixPtr[1][1] - matrixPtr[1][0] * matrixPtr[0][1] );
    
    }
    matrix_inverse *= det_inv;
  }
  else
  {

    // use numerical recipes Inverse of a Matrix

    int i, j;
    double d;
    std::vector<double> indx(numRows);
    std::vector<double> col(numRows);
    CubitMatrix save_matrix = *this;

    CubitStatus rv = ludcmp(&indx[0], d);
    if(rv != CUBIT_SUCCESS)
      throw std::invalid_argument ("rv must equal CUBIT_SUCCESS");
    //assert(rv == CUBIT_SUCCESS);
    for (j=0; j<numRows; j++)
    {
      for(i=0; i<numRows; i++) 
      {
        col[i] = 0.0;
      }
      col[j] = 1.0;
      rv = lubksb(&indx[0], &col[0]);
      if(rv != CUBIT_SUCCESS)
      throw std::invalid_argument ("rv must equal CUBIT_SUCCESS");
      //assert(rv == CUBIT_SUCCESS);
      for (i=0; i<numRows; i++) 
      {
        matrix_inverse.set(i,j,col[i]);
      }
    }
    *this = save_matrix;
  }
  
  return matrix_inverse;
}
bool CubitMatrix::is_equal ( const CubitMatrix &  other,
double  tol = 1e-8 
) const

Definition at line 1207 of file CubitMatrix.cpp.

{
  if ( numRows != other.numRows ) return false;
  if ( numCols != other.numCols ) return false;

  for (int i=0; i<numRows; i++)
  {
    for (int j=0; j<numCols; j++)
    {
      double diff = fabs(matrixPtr[i][j] - other.matrixPtr[i][j]);
      if(diff > tol)
        return false;
    }
  }
  return true;
}
bool CubitMatrix::is_identity ( double  tol = 1e-8) const

Definition at line 1183 of file CubitMatrix.cpp.

{
  bool ident = true;

  for (int i=0; i<numRows && ident; i++)
  {
    for (int j=0; j<numCols && ident; j++)
    {
      if (i == j)
      {
        if( fabs(matrixPtr[i][j] - 1.0) > tol)
          ident = false;
      }
      else
      {
        if(matrixPtr[i][j] > tol)
          ident = false;
      }
    }
  }

  return ident;
}
CubitStatus CubitMatrix::lubksb ( double *  indx,
double *  b 
)

Definition at line 1152 of file CubitMatrix.cpp.

{
  int i, j, ii, ip;
  double sum;

  // do the forward substitution

  ii = -1;
  for (i=0; i<numRows; ++i){
    ip = (int)indx[i];
    sum = b[ip];
    b[ip] = b[i];
    if (ii >= 0)
      for (j=ii; j<=i-1; ++j)
        sum -= matrixPtr[i][j]*b[j];
    else if (sum)
      ii = i;
    b[i] = sum;
  }

  // do the back substitution

  for (i=numRows-1; i>=0; --i){
    sum = b[i];
    for (j=i+1; j<numRows; ++j)
      sum -= matrixPtr[i][j]*b[j];
    b[i] = sum/matrixPtr[i][i];  // store a component of solution
  }
  return CUBIT_SUCCESS;
}
CubitStatus CubitMatrix::ludcmp ( double *  indx,
double &  d 
)

Definition at line 1078 of file CubitMatrix.cpp.

{
  int i, j, k, imax = -1;
  double big, tmp, sum;
  double *vv = new double [numRows];
  if (!vv) {
    return CUBIT_FAILURE;
  }

  d = 1.0; // no row interchanges yet

  // loop over rows to get implicit scale info

  for (i=0; i<numRows; ++i){
    big = 0.0;
    for (j=0; j<numRows; ++j)
      if ((tmp = fabs(matrixPtr[i][j])) > big)
        big = tmp;
    if (big == 0.0) {
        // note (vvyas, 3/2006): corrected array deletion
        // delete vv;
      delete [] vv;
      return CUBIT_FAILURE;
    }
    vv[i] = 1.0/big;
  }

  // loop over columns-Crout's method

  for (j=0; j<numRows; ++j){
    for (i=0; i<j; ++i){
      sum = matrixPtr[i][j];
      for (k=0; k<i; ++k)
        sum -= matrixPtr[i][k]*matrixPtr[k][j];
      matrixPtr[i][j] = sum;
    }
    big = 0.0;                         // initialize pivot search
    for (i=j; i<numRows; ++i){
      sum = matrixPtr[i][j];
      for (k=0; k<j; ++k)
        sum -= matrixPtr[i][k]*matrixPtr[k][j];
      matrixPtr[i][j] = sum;
      if ((tmp = vv[i]*fabs(sum)) > big) {
        big = tmp;
        imax = i;
      }
    }
    if (j != imax) {                   // do we need to change rows
      for (k=0; k<numRows; ++k) {
        tmp = matrixPtr[imax][k];
        matrixPtr[imax][k] = matrixPtr[j][k];
        matrixPtr[j][k] = tmp;
      }
      d = -d;
      vv[imax] = vv[j];
    }
    indx[j] = imax;
    if (matrixPtr[j][j] == 0.0) matrixPtr[0][0] = 1.0e-20;
    if (j != numRows-1) {             // divide by the pivot element
      tmp = 1.0/matrixPtr[j][j];
      for (i=j+1; i<numRows; ++i)
        matrixPtr[i][j] *= tmp;
    }
  }                                   // go back for next column

    // note (vvyas 3/2006): corrected array deletion
    // delete vv;
  delete [] vv;
  
  return CUBIT_SUCCESS;
}
int CubitMatrix::num_cols ( ) const [inline]

Definition at line 66 of file CubitMatrix.hpp.

{ return numCols; }
int CubitMatrix::num_rows ( ) const [inline]

Definition at line 65 of file CubitMatrix.hpp.

{ return numRows; }
CubitMatrix CubitMatrix::operator* ( const CubitMatrix &  matrix) const

Reimplemented in CubitTransformMatrix.

Definition at line 217 of file CubitMatrix.cpp.

{
    // Check that we can multiply them.
  if(numCols != matrix.num_rows())
    throw std::invalid_argument ("# of columns in first MUST match # of rows of second");
  //assert( numCols == matrix.num_rows() );
  
  CubitMatrix return_matrix( numRows, matrix.num_cols() );
  
  for( int ii = 0; ii < numRows; ii++ )
  {
    for( int jj = 0; jj < matrix.num_cols(); jj++ )
    {
      double temp = 0.0;
      for( int kk = 0; kk < numCols; kk++ )
      {
        //temp += get( ii, kk ) * matrix.get( kk, jj );
        temp += matrixPtr[ii][kk] * matrix.matrixPtr[kk][jj];
      }
      return_matrix.set( ii, jj, temp );
    }
  }
  return return_matrix;
}
CubitVector CubitMatrix::operator* ( const CubitVector &  vector) const

Reimplemented in CubitTransformMatrix.

Definition at line 245 of file CubitMatrix.cpp.

{
  // Check that we can multiply them.
  if(numCols!=3)
  { 
    throw std::invalid_argument("Matrix must have 3 columns");
  }
  //assert( numCols == 3 );
  
  double vec1[3];
  double vec2[3];
  
  vec2[0] = vector.x();
  vec2[1] = vector.y();
  vec2[2] = vector.z();
  
  for( int row = 0; row < numRows; row++ )
  {
    vec1[row] = 0.0;
    for( int col = 0; col < numCols; col++ )
    {
      vec1[row] += ( matrixPtr[row][col] * vec2[col] );
    }
  }
  
  return CubitVector( vec1[0],  vec1[1], vec1[2] );
}
std::vector< double > CubitMatrix::operator* ( const std::vector< double > &  vector) const

Definition at line 274 of file CubitMatrix.cpp.

{
    // Check that we can multiply them.
  if(numCols != (int) vector.size())
    throw std::invalid_argument ("Columns of Matrix do not match vector size");
  //assert( numCols == vector.size() );
  
  std::vector<double> return_vec( numRows );
  
  for( int row = 0; row < numRows; row++ )
  {
    return_vec[row] = 0.0;
    for( int col = 0; col < numCols; col++ )
    {
      return_vec[row] += ( matrixPtr[row][col] * vector[col] );
    }
  }
  
  return return_vec;
}
CubitMatrix CubitMatrix::operator* ( double  val) const

Reimplemented in CubitTransformMatrix.

Definition at line 296 of file CubitMatrix.cpp.

{
  CubitMatrix matrix( numRows, numCols );
  
  for( int row = 0; row < numRows; row++ )
  {
    for( int col = 0; col < numCols; col++ )
    {
      matrix.set( row, col,( matrixPtr[row][col] * val ) );
    }
  }
  return matrix;
}
CubitMatrix & CubitMatrix::operator*= ( const double  multiplier)

Definition at line 378 of file CubitMatrix.cpp.

{
  for( int ii = 0; ii < numRows; ii++ )
  {
    for( int jj = 0; jj < numCols; jj++ )
    {
      matrixPtr[ii][jj] *= multiplier;
    }
  }
  return *this;
}
CubitMatrix CubitMatrix::operator+ ( const CubitMatrix &  matrix) const

Definition at line 349 of file CubitMatrix.cpp.

{
   CubitMatrix return_matrix( numRows, numCols );
   
   for( int ii = 0; ii < numRows; ii++ )
   {
      for( int jj = 0; jj < numCols; jj++ )
      {
         return_matrix.set( ii, jj, matrixPtr[ii][jj] +
                            matrix.get( ii, jj ));
      }
   }

   return return_matrix;
}
CubitMatrix & CubitMatrix::operator+= ( const CubitMatrix &  matrix)

Definition at line 366 of file CubitMatrix.cpp.

{
  for( int ii = 0; ii < numRows; ii++ )
  {
    for( int jj = 0; jj < numCols; jj++ )
    {
      matrixPtr[ii][jj] += matrix.get( ii, jj );
    }
  }
  return *this;
}
CubitMatrix CubitMatrix::operator- ( const CubitMatrix &  matrix) const

Definition at line 331 of file CubitMatrix.cpp.

{
  CubitMatrix return_matrix( numRows, numCols );
  
  for( int ii = 0; ii < numRows; ii++ )
  {
    for( int jj = 0; jj < numCols; jj++ )
    {
      return_matrix.set( ii, jj, matrixPtr[ii][jj] -
                         matrix.get( ii, jj ));
    }
  }
  
  return return_matrix;
}
CubitMatrix CubitMatrix::operator/ ( double  val) const

Definition at line 311 of file CubitMatrix.cpp.

{
  if(val==0)
    throw std::invalid_argument("Cannot Divide by Zero");
  //assert( val != 0 );
  CubitMatrix matrix( numRows, numCols );
  
  for( int ii = 0; ii < numRows; ii++ )
  {
    for( int jj = 0; jj < numCols; jj++ )
    {
      matrix.set( ii, jj,( matrixPtr[ii][jj] / val ) );
    }
  }
  return matrix;
}
CubitMatrix CubitMatrix::operator= ( const CubitMatrix &  matrix)

Definition at line 189 of file CubitMatrix.cpp.

{
  int i, j;
  
  if (numRows != matrix.num_rows() ||
      numCols != matrix.num_cols())
  {
    delete [] matrixPtr;
    delete [] matrixMem;
    
    numRows = matrix.num_rows();
    numCols = matrix.num_cols();
    matrixPtr = new double*[numRows];
    matrixMem = new double[numRows*numCols];
    for (i = 0; i < numRows; i++)
      matrixPtr[i] = &matrixMem[i*numCols];
  }
  
  for(i = 0; i < numRows; i++ )
    for(j = 0; j < numCols; j++)
      matrixPtr[i][j] = matrix.get(i, j);
  
  return *this;
}
void CubitMatrix::plus_identity ( )

Definition at line 1224 of file CubitMatrix.cpp.

{
  for (int i=0; i<numRows; i++)
  {
    if ( i == numCols ) break;
    matrixPtr[i][i] += 1.0;
  }
}
CubitBoolean CubitMatrix::positive_definite ( ) const

Definition at line 514 of file CubitMatrix.cpp.

{

  if ( matrixPtr[0][0] <= 0. ) { return CUBIT_FALSE; }

  double det2x2 = matrixPtr[0][0] * matrixPtr[1][1] - matrixPtr[1][0] * matrixPtr[0][1];
  if ( det2x2 <= 0. ) { return CUBIT_FALSE; }

  if ( determinant() <= 0. ) { return CUBIT_FALSE; }
  
  return CUBIT_TRUE;
}
void CubitMatrix::print_matrix ( ) const

Definition at line 159 of file CubitMatrix.cpp.

{
  printf( "\n\n" );
  for( int row = 0; row < numRows; row++ )
  {
    for( int col = 0; col < numCols; col++ )
       PRINT_INFO("%25.15f", matrixPtr[row][col]);
    PRINT_INFO("\n");
  }
}
void CubitMatrix::print_matrix ( char *  filename) const

Definition at line 169 of file CubitMatrix.cpp.

{
  CubitFile fp( filename, "w" );
  if ( !fp )
  {
    printf( "CubitMatrix::print_matrix - Unable to open %s for writing\n",
            filename );
    return;
  }

  for( int row = 0; row < numRows; row++ )
  {
    for( int col = 0; col < numCols; col++ )
       fprintf( fp.file(), "%20.15f", matrixPtr[row][col] );
    fprintf( fp.file(), "\n" );
  }
}
int CubitMatrix::rank ( void  ) const

Definition at line 766 of file CubitMatrix.cpp.

{
  const double tol = 1E-12;
  int rank = 0;
  CubitMatrix tmp = *this;

  int irow;
  for ( irow = 0; irow < numRows; irow++ )
  {
    // make sure tmp[irow][irow] is non-zero.  If it isn't, swap a row to
    // make it so.
    double val = tmp.get(irow,irow);
    if ( fabs(val) < tol )
    {
      bool found = false;
      for ( int i = irow+1; i < numRows; i++ )
      {
        if ( fabs(tmp.get(i,irow)) > 1E-4 )
        {
          // swap row (irow) with row (irow+i).
          for ( int icol = 0; icol < numCols; icol++ )
          {
            double tmp1 = tmp.get(irow, icol);
            double tmp2 = tmp.get(i, icol);
            tmp.set(irow, icol, tmp2 );
            tmp.set(i, icol, tmp1 );
            found = true;
          }
        }
        if ( found ) break;
      }
      val = tmp.get(irow,irow);
    }
    if ( fabs(val) < tol )
      continue;

    rank++;

    for ( int icol = 0; icol < numCols; icol++ )
    {
      double col_val = tmp.get(irow, icol);
      tmp.set(irow,icol, col_val/val );
    }

    for ( int jrow = irow+1; jrow < numRows; jrow++ )
    {
      val = tmp.get(jrow,irow);
      if ( fabs(val) < tol )
        continue;

      for ( int icol = 0; icol < numCols; icol++ )
      {
        double tmp1 = tmp.get(jrow,icol) / val;
        tmp1 -= tmp.get(irow, icol);
        tmp.set(jrow,icol,tmp1 );
      }
    }
  }
  return rank;
}
void CubitMatrix::reset_size ( const int  n,
const int  m,
double  default_value = 0.0 
)

Definition at line 136 of file CubitMatrix.cpp.

{
  if ( matrixPtr ) delete [] matrixPtr;
  if ( matrixMem ) delete [] matrixMem;

  numRows = n;
  numCols = m;
  
  matrixMem = new double[numRows*numCols];
  matrixPtr = new double *[n];
  
  int ii;
  for(  ii = 0; ii < n; ii++ )
  {
    matrixPtr[ii] = &matrixMem[ii*numCols];
  }
  
    // Initialize matrix to zeros.     
  for( ii = 0; ii < n; ii++ )
     for( int jj = 0 ; jj < m; jj++ )
        matrixPtr[ii][jj] = default_value;
}
void CubitMatrix::set ( const int  n,
const int  m,
const double  val 
) [inline]

Definition at line 81 of file CubitMatrix.hpp.

    {
      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;
    }
void CubitMatrix::set_to_identity ( )

Definition at line 392 of file CubitMatrix.cpp.

{
  for (int i = numRows; i--; )
    for (int j = numCols; j--; )
    {
      if (i == j)
        matrixPtr[i][j] = 1;
      else
        matrixPtr[i][j] = 0;
    }
}
void CubitMatrix::solve ( CubitVector &  b,
const CubitVector &  pivot 
)

Definition at line 903 of file CubitMatrix.cpp.

{
  double rhs[3];
  rhs[0] = b.x();
  rhs[1] = b.y();
  rhs[2] = b.z();

  double pvt[3];
  pvt[0] = pivot.x();
  pvt[1] = pivot.y();
  pvt[2] = pivot.z();

  int j;
  const int n=3;
  for ( int k=0; k<n-1; k++ ) 
  {
     j=(int)pvt[k];
     if ( j != k ) 
     {
        double swap = rhs[k];
        rhs[k] = rhs[j];
        rhs[j] = swap;
     }

     for ( int i=k+1; i<n; i++ ) 
     {
        rhs[i] -= matrixPtr[i][k] * rhs[k];
     }

  }

  rhs[n-1] /= matrixPtr[n-1][n-1];

  for ( int i=n-2; i>-1; i-- ) 
  {
     double sum=0.;
     for ( j=i+1; j<n; j++ ) 
     {
        sum += matrixPtr[i][j] * rhs[j];
     }
     rhs[i] = ( rhs[i] - sum ) / matrixPtr[i][i];
  }

  b.set( rhs[0], rhs[1], rhs[2] );

}
CubitStatus CubitMatrix::solveNxN ( CubitMatrix &  rhs,
CubitMatrix &  coef 
)

Definition at line 1013 of file CubitMatrix.cpp.

{
  if (numRows != rhs.num_rows() ||
    numRows != numCols) {
    return CUBIT_FAILURE;
  }
  int i,j;
  double d;
  double *indx = new double [numRows];
  double *b = new double [numRows];

  CubitStatus status = ludcmp(indx, d);
  if (status == CUBIT_SUCCESS)
  {
    coef.reset_size( rhs.num_rows(), rhs.num_cols(), 0.0 );
    for ( j = 0; j < rhs.num_cols(); j++ )
    {
      for(i=0; i<numRows; i++)
      {
        b[i] = rhs.get(i,j);
      }
      status = lubksb(indx, b);
      for (i=0; i<numRows; i++)
      {
        coef.set(i,j,b[i]);
      }
    }
  }
  delete [] indx;
  delete [] b;
  return status;
}
CubitStatus CubitMatrix::solveNxN ( const std::vector< double > &  rhs,
std::vector< double > &  coef 
)

Definition at line 1046 of file CubitMatrix.cpp.

{
  if (numRows != (int) rhs.size() ||
    numRows != numCols) {
    return CUBIT_FAILURE;
  }
  int i;
  double d;
  double *indx = new double [numRows];
  double *b = new double [numRows];

  CubitStatus status = ludcmp(indx, d);
  if (status == CUBIT_SUCCESS)
  {
    coef.clear();
    for(i=0; i<numRows; i++)
    {
      b[i] = rhs[i];
    }
    status = lubksb(indx, b);
    for (i=0; i<numRows; i++)
    {
      coef.push_back( b[i] );
    }
  }
  delete [] indx;
  delete [] b;
  return status;
}
CubitMatrix CubitMatrix::sub_matrix ( const int  row,
const int  col 
) const

Definition at line 597 of file CubitMatrix.cpp.

{
  CubitMatrix matrix (numRows - 1, numCols - 1);
  
  int copy_row = 0;
  for (int source_row = 0; source_row < numRows; source_row++)
  {
    if (source_row != row)
    {
      int copy_col = 0;
      for (int source_col = 0; source_col < numCols; source_col++)
      {
        if (source_col != col)
        {
          matrix.set (copy_row, copy_col, matrixPtr[source_row][source_col]);
          copy_col++;
        }
      }
      copy_row++;
    }
  }
  
  return matrix;
}
void CubitMatrix::sub_matrix ( const std::vector< bool > &  rows_to_include,
const std::vector< bool > &  cols_to_include,
CubitMatrix &  submatrix 
)

Definition at line 625 of file CubitMatrix.cpp.

{
  if(numRows != (int) rows_to_include.size())
    throw std::invalid_argument ("rows_to_include size must match numRows");
  //assert( numRows == rows_to_include.size() );
  if(numCols != (int) cols_to_include.size())
    throw std::invalid_argument ("cols_to_include size must match numCols");
  //assert( numCols == cols_to_include.size() );

  int i;
  int nrow = 0, ncol = 0;
  for ( i = 0; i < numRows; i++ ) if ( rows_to_include[i] ) nrow++;
  for ( i = 0; i < numCols; i++ ) if ( cols_to_include[i] ) ncol++;
  submatrix.reset_size( nrow, ncol, 0.0 );

  for ( int r = 0, new_r = 0; r < numRows; r++ )
  {
    if ( !rows_to_include[r] ) continue;
    for ( int c = 0, new_c = 0; c < numCols; c++ )
    {
      if ( !cols_to_include[c] ) continue;
      submatrix.set( new_r, new_c, get(r,c) );
      new_c++;
    }
    new_r++;
  }
}
CubitMatrix CubitMatrix::symm ( ) const
CubitMatrix CubitMatrix::transpose ( ) const

Definition at line 580 of file CubitMatrix.cpp.

{
  CubitMatrix return_matrix( numCols, numRows );
  
  for( int ii = 0; ii < numRows; ii++ )
  {
    for( int jj = 0; jj < numCols; jj++ )
    {
      return_matrix.set( jj, ii, matrixPtr[ii][jj] );
    }
  }
  
  return return_matrix;
}
bool CubitMatrix::vector_outer_product ( const CubitVector &  vec1,
const CubitVector &  vec2 
)

Definition at line 1233 of file CubitMatrix.cpp.

{             
  if ( numRows != 3 || numCols != 3 )
    return false;

    // Initialize the matrix elements using otimes (outer product)
  matrixPtr[0][0] = vec1.x() * vec2.x();
  matrixPtr[1][0] = vec1.y() * vec2.x();
  matrixPtr[2][0] = vec1.z() * vec2.x();
  matrixPtr[0][1] = vec1.x() * vec2.y();
  matrixPtr[1][1] = vec1.y() * vec2.y();
  matrixPtr[2][1] = vec1.z() * vec2.y();
  matrixPtr[0][2] = vec1.x() * vec2.z();
  matrixPtr[1][2] = vec1.y() * vec2.z();
  matrixPtr[2][2] = vec1.z() * vec2.z();
  return true;
}

Member Data Documentation

double* CubitMatrix::matrixMem [private]

Definition at line 171 of file CubitMatrix.hpp.

double** CubitMatrix::matrixPtr [private]

Definition at line 170 of file CubitMatrix.hpp.

int CubitMatrix::numCols [private]

Definition at line 173 of file CubitMatrix.hpp.

int CubitMatrix::numRows [private]

Definition at line 172 of file CubitMatrix.hpp.

The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines