matrix.c

开放源码的编译器open watcom 1.6.0版的源代码

  int i, j;
  for (i = 0; i < this->num_rows; i++)
    for (j = 0; j < this->num_cols; j++)
      this->data[i][j] -= m.data[i][j];
  return *this;
}

// operator* -- Non Destructive matrix multiply 
//               num_cols of first matrix must match num_rows of second matrix.
// Input:        two matrix references
// Output:       New matrix containing the product.

template<class Type>
CoolEnvelope< CoolMatrix<Type> > operator* (const CoolMatrix<Type>& m1,
                             const CoolMatrix<Type>& m2) {
  if (m1.num_cols != m2.num_rows)               // dimensions do not match?
    m1.dimension_error ("operator*=", "Type", 
                        m1.num_rows, m1.num_cols, m2.num_rows, m2.num_cols);
  CoolMatrix<Type> temp(m1.num_rows, m2.num_cols); // Temporary to store product
  for (int i = 0; i < m1.num_rows; i++) {       // For each row
    for (int j = 0; j < m2.num_cols; j++) {     // For each element in column
      Type sum = 0;
      for (int k = 0; k < m1.num_cols; k++)     // Loop over column values
        sum += (m1.data[i][k] * m2.data[k][j]); // Multiply
      temp(i,j) = sum;
    }
  }
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}

// operator- -- Non-destructive matrix negation. a = -b;
// Input:       this*
// Output:      New matrix 

template<class Type> 
CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::operator- () const {
  CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  int i, j;
  for (i = 0; i < this->num_rows; i++)
    for (j = 0; j < this->num_cols; j++)
      temp.data[i][j] = - this->data[i][j];
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}

// operator+ -- Non-destructive matrix addition of a scalar.
// Input:       this*, scalar value
// Output:      New matrix 

template<class Type> 
CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::operator+ (const Type& value) const {
  CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  for (int i = 0; i < this->num_rows; i++)      // For each row
    for (int j = 0; j < this->num_cols; j++)    // For each element in column
      temp.data[i][j] = (this->data[i][j] + value); // Add scalar
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}


// operator* -- Non-destructive matrix multiply by a scalar.
// Input:       this*, scalar value
// Output:      New matrix 

template<class Type> 
CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::operator* (const Type& value) const {
  CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  for (int i = 0; i < this->num_rows; i++)      // For each row
    for (int j = 0; j < this->num_cols; j++)    // For each element in column
      temp.data[i][j] = (this->data[i][j] * value); // Multiply
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}


// operator/ -- Non-destructive matrix division by a scalar.
// Input:       this*, scalar value
// Output:      New matrix 

template<class Type> 
CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::operator/ (const Type& value) const {
  CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  for (int i = 0; i < this->num_rows; i++)      // For each row
    for (int j = 0; j < this->num_cols; j++)    // For each element in column
      temp.data[i][j] = (this->data[i][j] / value); // Divide
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}


////--------------------------- Additions ------------------------------------


// transpose -- Return the transpose of this matrix.
// Input:       this*
// Ouput:       New matrix

template<class Type> 
CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::transpose() const {
  CoolMatrix<Type> temp(this->num_cols, this->num_rows);
  int i, j;
  for (i = 0; i < this->num_cols; i++)
    for (j = 0; j < this->num_rows; j++)
      temp.data[i][j] = this->data[j][i];
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}


// abs -- Return the matrix of the absolute values.
// Input:       this*
// Ouput:       New matrix

template<class Type> 
CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::abs() const {
  CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  int i, j;
  for (i = 0; i < this->num_rows; i++)
    for (j = 0; j < this->num_cols; j++)
      if (this->data[i][j] < 0)
        temp.data[i][j] = - this->data[i][j];
      else
        temp.data[i][j] = this->data[i][j];
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}

// sign -- Return the matrix whose elements are either -1,1 or 0
// depending on whether the corresponding values are negative, positive, or 0.
// Input:       this*
// Ouput:       New matrix

template<class Type> 
CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::sign() const {
  CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  int i, j;
  for (i = 0; i < this->num_rows; i++)
    for (j = 0; j < this->num_cols; j++)
      if (this->data[i][j] == 0)                        // test fuzz equality to 0
        temp.data[i][j] = 0;                    // first.
      else
        if (this->data[i][j] < 0)
          temp.data[i][j] = -1;
        else
          temp.data[i][j] = 1;
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}

// element_product -- return the matrix whose elements are the products 
// Input:     2 matrices m1, m2 by reference
// Output:    New matrix, whose elements are m1[ij]*m2[ij].

template<class Type>
CoolEnvelope< CoolMatrix<Type> > element_product (const CoolMatrix<Type>& m1, const CoolMatrix<Type>& m2) {
  if (m1.num_rows != m2.num_rows || m1.num_cols != m2.num_cols) // Size?
    m1.dimension_error ("element_product", "Type", 
                        m1.num_rows, m1.num_cols, m2.num_rows, m2.num_cols);
  CoolMatrix<Type> temp(m1.num_rows, m1.num_cols);
  int i, j;
  for (i = 0; i < m1.num_rows; i++)
    for (j = 0; j < m1.num_cols; j++)
      temp.data[i][j] = m1.data[i][j] * m2.data[i][j];
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}

// element_quotient -- return the matrix whose elements are the quotients 
// Input:     2 matrices m1, m2 by reference
// Output:    New matrix, whose elements are m1[ij]/m2[ij].

template<class Type>
CoolEnvelope< CoolMatrix<Type> > element_quotient (const CoolMatrix<Type>& m1, const CoolMatrix<Type>& m2) {
  if (m1.num_rows != m2.num_rows || m1.num_cols != m2.num_cols) // Size?
    m1.dimension_error ("element_quotient", "Type", 
                        m1.num_rows, m1.num_cols, m2.num_rows, m2.num_cols);
  CoolMatrix<Type> temp(m1.num_rows, m1.num_cols);
  int i, j;
  for (i = 0; i < m1.num_rows; i++)
    for (j = 0; j < m1.num_cols; j++)
      temp.data[i][j] = m1.data[i][j] / m2.data[i][j];
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}

// update -- replace a submatrix of this, by the actual argument.
// Input:       *this, starting corner specified by top and left.
// Ouput:       mutated reference.

template<class Type> 
CoolMatrix<Type>& CoolMatrix<Type>::update (const CoolMatrix<Type>& m, 
                                            unsigned int top, unsigned int left) {
  unsigned int bottom = top + m.num_rows;
  unsigned int right = left + m.num_cols;
  if (this->num_rows < bottom || this->num_cols < right)
    this->dimension_error ("update", "Type", 
                           bottom-top, right-left, m.num_rows, m.num_cols);
  int i, j;
  for (i = top; i < bottom; i++)
    for (j = left; j < right; j++)
      this->data[i][j] = m.data[i-top][j-left];
  return *this;
}


// extract -- Return a submatrix specified by the top-left corner and size.
// Input:       *this, starting corner specified by top and left, and size.
// Ouput:       new matrix

template<class Type> 
CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::extract (unsigned int rows, unsigned int cols, unsigned int top, unsigned int left) const{
  unsigned int bottom = top + rows;
  unsigned int right = left + cols;
  if ((this->num_rows < bottom) || (this->num_cols < right))
    this->dimension_error ("extract", "Type", 
                           bottom-top, right-left, rows, cols);
  CoolMatrix<Type> temp(rows, cols);
  for (int i = 0; i < rows; i++)                // actual copy of all elements
    for (int j = 0; j < cols ; j++)             // in submatrix
      temp.data[i][j] = data[top+i][left+j];
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}

// determinant -- Determinant of a square matrix using Cramer's rule.
//              Signal Error exception if the matrix is not square.

template<class Type>
Type CoolMatrix<Type>::determinant () const {
  if (this->num_rows != this->num_cols || this->num_rows < 2)
    this->dimension_error ("determinant", "Type",
                           this->num_rows, this->num_cols, 
                           this->num_rows, this->num_cols);
  int n = this->num_rows, r, i, j;
  Type det = 0, prod;
  if (n == 2) {
    det = (this->data[0][0] * this->data[1][1] - // border case of 2x2 matrix
           this->data[0][1] * this->data[1][0]);
  } else {
    for (r = 0; r < n; r++) {                   // compute sum of products
      prod = 1;                                 // along diagonals
      for (i = r, j = 0; i < n; i++, j++)       // top-left to bottom-right
        prod *= this->data[i][j];
      for (i = 0; j < n; i++, j++)
        prod *= this->data[i][j];
      det += prod;                              // coeft = +1
    }
    int e = n-1;                                // index of last row/col
    for (r = 0; r < n; r++) {                   // compute sum of products
      prod = 1;                                 // along diagonals
      for (i = r, j = e; i < n; i++, j--)       // top-right to bottom-left
        prod *= this->data[i][j];
      for (i=0; j >= 0; i++, j--)
        prod *= this->data[i][j];
    det -= prod;                                // coeft = -1
    }
  }
  return det;
}

// dot_product -- Return the dot product of the row or column vectors
// Input:       2 vectors by reference
// Ouput:       dot product value

template<class Type>
  Type dot_product (const CoolMatrix<Type>& v1, const CoolMatrix<Type>& v2) {
    if (v1.num_rows != v2.num_rows || v1.num_cols != v2.num_cols) // Size?
      v1.dimension_error ("dot_product", "Type", 
                          v1.num_rows, v1.num_cols, v2.num_rows, v2.num_cols);
    Type dot = 0;
    int i, j;
    for (i = 0; i < v1.num_rows; i++)
      for (j = 0; j < v1.num_cols; j++)         // generalized dot-product
        dot += v1.data[i][j] * v2.data[i][j];   // of matrices
    return dot;
}

// cross_2d -- Return the 2X1 cross-product of 2 2d-vectors
// Input:       2 vectors by reference
// Ouput:       cross product value

template<class Type>
Type cross_2d (const CoolMatrix<Type>& v1, const CoolMatrix<Type>& v2) {
  if (v1.num_rows != v2.num_rows || v1.num_cols != v2.num_cols)
    v1.dimension_error ("cross_2d", "Type", 
                        v1.num_rows, v1.num_cols, v2.num_rows, v2.num_cols);
  CoolMatrix<Type>& m1 = (CoolMatrix<Type>&) v1; // cast away const.
  CoolMatrix<Type>& m2 = (CoolMatrix<Type>&) v2; 
  return (m1.x() * m2.y()                       // work for both col/row
          -                                     // representation.
          m1.y() * m2.x());
}

// cross_3d -- Return the 3X1 cross-product of 2 3d-vectors
// Input:       2 vectors by reference
// Ouput:       3d cross product vector

template<class Type>
CoolEnvelope< CoolMatrix<Type> > cross_3d (const CoolMatrix<Type>& v1, const CoolMatrix<Type>& v2) {
  if (v1.num_rows != v2.num_rows || v1.num_cols != v2.num_cols)
    v1.dimension_error ("cross_3d", "Type", 
                        v1.num_rows, v1.num_cols, v2.num_rows, v2.num_cols);
  CoolMatrix<Type> temp(v1.num_rows, v1.num_cols);
  CoolMatrix<Type>& m1 = (CoolMatrix<Type>&) v1; // cast away const.
  CoolMatrix<Type>& m2 = (CoolMatrix<Type>&) v2; 
  temp.x() = m1.y() * m2.z() - m1.z() * m2.y(); // work for both col/row
  temp.y() = m1.z() * m2.x() - m1.x() * m2.z(); // representation
  temp.z() = m1.x() * m2.y() - m1.y() * m2.x();
  CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  return result;                                         // copy of envelope
}

//## hack to workaround BC++ 3.1 Envelope bug
#undef CoolEnvelope