Line data    Source code

       1             : /**
       2             :  * @file
       3             :  * \brief Defines a Matrix class
       4             :  *
       5             :  * @author Alicia Klinvex
       6             :  */
       7             : 
       8             : #include "Morpheus_Matrix.h"
       9             : #include <cassert>
      10             : #include <cmath>
      11             : #include <iostream>
      12             : 
      13             : namespace Morpheus {
      14             : 
      15           3 : Matrix::Matrix(const int nrows, const int ncols)
      16             : {
      17           3 :   nrows_ = nrows;
      18           3 :   ncols_ = ncols;
      19             : 
      20             :   // Make sure the dimensions make sense
      21           3 :   assert(nrows_ > 0);
      22           3 :   assert(ncols_ > 0);
      23             : 
      24             :   // Allocate memory for the data
      25           3 :   data_ = new double*[nrows_];
      26          18 :   for(int r=0; r<nrows_; r++)
      27          15 :     data_[r] = new double[ncols_];
      28           3 : }
      29             : 
      30             : 
      31           3 : Matrix::~Matrix()
      32             : {
      33             :   // Free all the memory we allocated
      34          18 :   for(int r=0; r<nrows_; r++)
      35          15 :     delete[] data_[r];
      36             : 
      37           3 :   delete[] data_;
      38           3 : }
      39             : 
      40             : 
      41         325 : double& Matrix::operator()(const int row, const int col)
      42             : {
      43         325 :   return data_[row][col];
      44             : }
      45             : 
      46             : 
      47           0 : void Matrix::multiply(const Vector& X, Vector& Y) const
      48             : {
      49             :   // Make sure the dimensions are consistent
      50           0 :   assert(X.getNumElements() == ncols_);
      51           0 :   assert(Y.getNumElements() == nrows_);
      52             : 
      53           0 :   for(int r=0; r<nrows_; r++)
      54             :   {
      55           0 :     Y[r] = 0;
      56           0 :     for(int c=0; c<ncols_; c++)
      57             :     {
      58           0 :       Y[r] = Y[r] + data_[r][c]*X[c];
      59             :     }
      60             :   }
      61           0 : }
      62             : 
      63             : 
      64           1 : void Matrix::multiply(const Matrix& X, Matrix& Y) const
      65             : {
      66             :   // Make sure the dimensions are consistent
      67           1 :   assert(nrows_ == Y.nrows_);
      68           1 :   assert(ncols_ == X.nrows_);
      69           1 :   assert(X.ncols_ == Y.ncols_);
      70             : 
      71           6 :   for(int r=0; r<Y.nrows_; r++)
      72             :   {
      73          30 :     for(int c=0; c<Y.ncols_; c++)
      74             :     {
      75          25 :       Y(r,c) = 0;
      76         150 :       for(int k=0; k<ncols_; k++)
      77             :       {
      78         125 :         Y(r,c) = Y(r,c) + data_[r][k] * X.data_[k][c];
      79             :       }
      80             :     }
      81             :   }
      82           1 : }
      83             : 
      84             : 
      85           1 : bool Matrix::isSymmetric() const
      86             : {
      87           1 :   if(nrows_ != ncols_)
      88           0 :     return false;
      89             : 
      90           6 :   for(int r=0; r<nrows_; r++)
      91             :   {
      92          30 :     for(int c=0; c<ncols_; c++)
      93             :     {
      94          25 :       if(data_[r][c] != data_[c][r])
      95           0 :         return false;
      96             :     }
      97             :   }
      98             : 
      99           1 :   return true;
     100             : }
     101             : 
     102             : 
     103           1 : bool Matrix::isUpperTriangular() const
     104             : {
     105           1 :   if(nrows_ != ncols_)
     106           0 :     return false;
     107             : 
     108           6 :   for(int r=0; r<nrows_; r++)
     109             :   {
     110          15 :     for(int c=r+1; c<ncols_; c++)
     111             :     {
     112          10 :       if(data_[r][c] != 0)
     113           0 :         return false;
     114             :     }
     115             :   }
     116             : 
     117           1 :   return true;
     118             : }
     119             : 
     120             : 
     121             : // Maximum absolute column sum
     122           0 : double Matrix::norm1() const
     123             : {
     124             :   double curColSum;
     125           0 :   double maxColSum = 0;
     126           0 :   for(int c=0; c<ncols_; c++)
     127             :   {
     128           0 :     curColSum = 0;
     129           0 :     for(int r=0; r<nrows_; r++)
     130             :     {
     131           0 :       curColSum = curColSum + data_[r][c];
     132             :     }
     133           0 :     if(curColSum > maxColSum)
     134           0 :       maxColSum = curColSum;
     135             :   }
     136           0 :   return maxColSum;
     137             : }
     138             : 
     139             : 
     140             : // Maximum absolute row sum
     141           0 : double Matrix::normInf() const
     142             : {
     143             :   double curRowSum;
     144           0 :   double maxRowSum = 0;
     145           0 :   for(int r=0; r<nrows_; r++)
     146             :   {
     147           0 :     curRowSum = 0;
     148           0 :     for(int c=0; c<ncols_; c++)
     149             :     {
     150           0 :       curRowSum = curRowSum + data_[r][c];
     151             :     }
     152           0 :     if(curRowSum > maxRowSum)
     153           0 :       maxRowSum = curRowSum;
     154             :   }
     155           0 :   return maxRowSum;
     156             : }
     157             : 
     158             : 
     159           0 : int Matrix::getNumRows() const
     160             : {
     161           0 :   return nrows_;
     162             : }
     163             : 
     164             : 
     165           0 : int Matrix::getNumCols() const
     166             : {
     167           0 :   return ncols_;
     168             : }
     169             : 
     170             : 
     171           0 : int Matrix::getNumEntries() const
     172             : {
     173           0 :   return (nrows_*ncols_);
     174             : }
     175             : 
     176             : 
     177           1 : bool Matrix::approxEqual(const Matrix& m, const double tol) const
     178             : {
     179           1 :   if(nrows_ != m.nrows_ || ncols_ != m.ncols_)
     180           0 :     return false;
     181             : 
     182           6 :   for(int r=0; r<nrows_; r++)
     183             :   {
     184          30 :     for(int c=0; c<ncols_; c++)
     185             :     {
     186          25 :       if(std::abs(data_[r][c] - m.data_[r][c]) > tol)
     187           0 :         return false;
     188             :     }
     189             :   }
     190           1 :   return true;
     191             : }
     192             : 
     193             : 
     194           0 : void Matrix::print() const
     195             : {
     196           0 :   std::cout << nrows_ << "x" << ncols_ << " Matrix\n";
     197           0 :   for(int r=0; r<nrows_; r++)
     198             :   {
     199           0 :     for(int c=0; c<ncols_; c++)
     200             :     {
     201           0 :       std::cout << data_[r][c] << " ";
     202             :     }
     203           0 :     std::cout << std::endl;
     204             :   }
     205           0 : }
     206             : 
     207           3 : } /* namespace Morpheus */