Logo Search packages:      
Sourcecode: octave2.9 version File versions

MatrixType.cc

/*

Copyright (C) 2006 David Bateman
Copyright (C) 2006 Andy Adler

Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.

Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, write to the
Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
Boston, MA 02110-1301, USA.

*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <vector>

#include "MatrixType.h"
#include "dMatrix.h"
#include "CMatrix.h"
#include "dSparse.h"
#include "CSparse.h"
#include "oct-spparms.h"

// FIXME There is a large code duplication here

MatrixType::MatrixType (void)
  : typ (MatrixType::Unknown), 
    sp_bandden (octave_sparse_params::get_bandden()),
    bandden (0), upper_band (0), 
    lower_band (0), dense (false), full (false), nperm (0), perm (0) { }

MatrixType::MatrixType (const MatrixType &a)
  : typ (a.typ), sp_bandden (a.sp_bandden), bandden (a.bandden), 
    upper_band (a.upper_band), lower_band (a.lower_band), 
    dense (a.dense), full (a.full), nperm (a.nperm)
{ 
  if (nperm != 0)
    {
      perm = new octave_idx_type [nperm];
      for (octave_idx_type i = 0; i < nperm; i++)
      perm[i] = a.perm[i];
    }
}

MatrixType::MatrixType (const Matrix &a)
  : typ (MatrixType::Unknown),
    sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
    dense (false), full (true), nperm (0), perm (0)
{
  octave_idx_type nrows = a.rows ();
  octave_idx_type ncols = a.cols ();
 
  if (ncols == nrows)
    {
      bool upper = true;
      bool lower = true;
      bool hermitian = true;

      for (octave_idx_type j = 0; j < ncols; j++)
      {
        if (j < nrows)
          {
            if (a.elem (j,j) == 0.)
            {
              upper = false;
              lower = false;
              hermitian = false;
              break;
            }
            if (a.elem (j,j) < 0.)
            hermitian = false;
          }      
        for (octave_idx_type i = 0; i < j; i++)
          if (lower && a.elem (i,j) != 0.)
            {
            lower = false;
            break;
            }
        for (octave_idx_type i = j+1; i < nrows; i++)
          {
            if (hermitian && a.elem (i, j) != a.elem (j, i))
            hermitian = false;
            if (upper && a.elem (i,j) != 0)
            upper = false;
          }
        if (!upper && !lower && !hermitian)
          break;
      }

      if (upper)
      typ = MatrixType::Upper;
      else if (lower)
      typ = MatrixType::Lower;
      else if (hermitian)
      typ = MatrixType::Hermitian;
      else if (ncols == nrows)
      typ = MatrixType::Full;
    }
  else
    typ = MatrixType::Rectangular;
}

MatrixType::MatrixType (const ComplexMatrix &a)
  : typ (MatrixType::Unknown),
    sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
    dense (false), full (true), nperm (0), perm (0)
{
  octave_idx_type nrows = a.rows ();
  octave_idx_type ncols = a.cols ();

  if (ncols == nrows)
    {
      bool upper = true;
      bool lower = true;
      bool hermitian = true;

      for (octave_idx_type j = 0; j < ncols; j++)
      {
        if (j < ncols)
          {
            if (imag(a.elem (j,j)) == 0. && 
              real(a.elem (j,j)) == 0.)
            {
              upper = false;
              lower = false;
              hermitian = false;
              break;
            }

            if (imag(a.elem (j,j)) != 0. || 
              real(a.elem (j,j)) < 0.)
                hermitian = false;
          }
        for (octave_idx_type i = 0; i < j; i++)
          if (lower && (real(a.elem (i,j)) != 0 || imag(a.elem (i,j)) != 0))
            {
            lower = false;
            break;
            }
        for (octave_idx_type i = j+1; i < nrows; i++)
          {
            if (hermitian && a.elem (i, j) != conj(a.elem (j, i)))
            hermitian = false;
            if (upper && (real(a.elem (i,j)) != 0 || 
                      imag(a.elem (i,j)) != 0))
            upper = false;
          }
        if (!upper && !lower && !hermitian)
          break;
      }

      if (upper)
      typ = MatrixType::Upper;
      else if (lower)
      typ = MatrixType::Lower;
      else if (hermitian)
      typ = MatrixType::Hermitian;
      else if (ncols == nrows)
      typ = MatrixType::Full;
    }
  else
    typ = MatrixType::Rectangular;
}

MatrixType::MatrixType (const SparseMatrix &a)
  : typ (MatrixType::Unknown),
    sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
    dense (false), full (false), nperm (0), perm (0)
{
  octave_idx_type nrows = a.rows ();
  octave_idx_type ncols = a.cols ();
  octave_idx_type nm = (ncols < nrows ? ncols : nrows);
  octave_idx_type nnz = a.nzmax ();

  if (octave_sparse_params::get_key ("spumoni") != 0.)
    (*current_liboctave_warning_handler) 
      ("Calculating Sparse Matrix Type");

  sp_bandden = octave_sparse_params::get_bandden();
  bool maybe_hermitian = false;
  typ = MatrixType::Full;

  if (nnz == nm)
    {
      matrix_type tmp_typ = MatrixType::Diagonal;
      octave_idx_type i;
      // Maybe the matrix is diagonal
      for (i = 0; i < nm; i++)
      {
        if (a.cidx(i+1) != a.cidx(i) + 1)
          {
            tmp_typ = MatrixType::Full;
            break;
          }
        if (a.ridx(i) != i)
          {
            tmp_typ = MatrixType::Permuted_Diagonal;
            break;
          }
      }
        
      if (tmp_typ == MatrixType::Permuted_Diagonal)
      {
        std::vector<bool> found (nrows);

        for (octave_idx_type j = 0; j < i; j++)
          found [j] = true;
        for (octave_idx_type j = i; j < nrows; j++)
          found [j] = false;
            
        for (octave_idx_type j = i; j < nm; j++)
          {
            if ((a.cidx(j+1) > a.cidx(j) + 1)  || 
              ((a.cidx(j+1) == a.cidx(j) + 1) && found [a.ridx(j)]))
            {
              tmp_typ = MatrixType::Full;
              break;
            }
            found [a.ridx(j)] = true;
          }
      }
      typ = tmp_typ;
    }

  if (typ == MatrixType::Full)
    {
      // Search for banded, upper and lower triangular matrices
      bool singular = false;
      upper_band = 0;
      lower_band = 0;
      for (octave_idx_type j = 0; j < ncols; j++)
      {
        bool zero_on_diagonal = false;
        if (j < nrows)
          {
            zero_on_diagonal = true;
            for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
            if (a.ridx(i) == j)
              {
                zero_on_diagonal = false;
                break;
              }
          }

        if (zero_on_diagonal)
          {
            singular = true;
            break;
          }

        if (a.cidx(j+1) != a.cidx(j))
          {
            octave_idx_type ru = a.ridx(a.cidx(j));
            octave_idx_type rl = a.ridx(a.cidx(j+1)-1);

            if (j - ru > upper_band)
            upper_band = j - ru;
              
            if (rl - j > lower_band)
            lower_band = rl - j;
          }
      }

      if (!singular)
      {
        bandden = double (nnz) /
          (double (ncols) * (double (lower_band) +
                         double (upper_band)) -
           0.5 * double (upper_band + 1) * double (upper_band) -
           0.5 * double (lower_band + 1) * double (lower_band));

        if (nrows == ncols && sp_bandden != 1. && bandden > sp_bandden)
          {
            if (upper_band == 1 && lower_band == 1)
            typ = MatrixType::Tridiagonal;
            else
            typ = MatrixType::Banded;

            octave_idx_type nnz_in_band = 
            (upper_band + lower_band + 1) * nrows -
            (1 + upper_band) * upper_band / 2 -
            (1 + lower_band) * lower_band / 2;
            if (nnz_in_band == nnz)
            dense = true;
            else 
            dense = false;
          }
        else if (upper_band == 0)
          typ = MatrixType::Lower;
        else if (lower_band == 0)
          typ = MatrixType::Upper;

        if (upper_band == lower_band && nrows == ncols)
          maybe_hermitian = true;
      }

      if (typ == MatrixType::Full)
      {
        // Search for a permuted triangular matrix, and test if
        // permutation is singular

        // FIXME
        // Perhaps this should be based on a dmperm algorithm
        bool found = false;

        nperm = ncols;
        perm = new octave_idx_type [ncols];

        for (octave_idx_type i = 0; i < ncols; i++)
          perm [i] = -1;

        for (octave_idx_type i = 0; i < nm; i++)
          {
            found = false;

            for (octave_idx_type j = 0; j < ncols; j++)
            {
              if ((a.cidx(j+1) - a.cidx(j)) > 0 && 
                  (a.ridx(a.cidx(j+1)-1) == i))
                {
                  perm [i] = j;
                  found = true;
                  break;
                }
            }

            if (!found)
            break;
          }

        if (found)
          {
            typ = MatrixType::Permuted_Upper;
            if (ncols > nrows)
            {
              octave_idx_type k = nrows;
              for (octave_idx_type i = 0; i < ncols; i++)
                if (perm [i] == -1)
                  perm[i] = k++;
            }
          }
        else if (a.cidx(nm) == a.cidx(ncols))
          {
            nperm = nrows;
            delete [] perm;
            perm = new octave_idx_type [nrows];
            OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, nrows);

            for (octave_idx_type i = 0; i < nrows; i++)
            {
              perm [i] = -1;
              tmp [i] = -1;
            }

            for (octave_idx_type j = 0; j < ncols; j++)
            for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
                perm [a.ridx(i)] = j;

            found = true;
            for (octave_idx_type i = 0; i < nm; i++)
            if (perm[i] == -1)
              {
                found = false;
                break;
              }
            else
              {
                tmp[perm[i]] = 1;
              }

            if (found)
            {
              octave_idx_type k = ncols;
              for (octave_idx_type i = 0; i < nrows; i++)
                {
                  if (tmp[i] == -1)
                  {
                    if (k < nrows)
                      {
                        perm[k++] = i;
                      }
                    else
                      {
                        found = false;
                        break;
                      }
                  }
                }
            }

            if (found)
            typ = MatrixType::Permuted_Lower;
            else
            {
              delete [] perm;
              nperm = 0;
            }
          }
        else
          {
            delete [] perm;
            nperm = 0;
          }
      }

      // FIXME
      // Disable lower under-determined and upper over-determined problems
      // as being detected, and force to treat as singular. As this seems
      // to cause issues
      if (((typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
         && nrows > ncols) ||
        ((typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
         && nrows < ncols))
      {
        typ = MatrixType::Rectangular;
        if (typ == MatrixType::Permuted_Upper ||
            typ == MatrixType::Permuted_Lower)
          delete [] perm;
        nperm = 0;
      }

      if (typ == MatrixType::Full && ncols != nrows)
      typ = MatrixType::Rectangular;

      if (maybe_hermitian && (typ == MatrixType::Full || 
                        typ == MatrixType::Tridiagonal || 
                        typ == MatrixType::Banded))
      {
        // Check for symmetry, with positive real diagonal, which
        // has a very good chance of being symmetric positive
        // definite..
        bool is_herm = true;

        for (octave_idx_type j = 0; j < ncols; j++)
          {
            bool diag_positive = false;

            for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
            {
              octave_idx_type ri = a.ridx(i);

              if (ri == j)
                {
                  if (a.data(i) == std::abs(a.data(i)))
                  diag_positive = true;
                  else
                  break;
                }
              else
                {
                  bool found = false;

                  for (octave_idx_type k = a.cidx(ri); k < a.cidx(ri+1); k++)
                  {
                    if (a.ridx(k) == j)
                      {
                        if (a.data(i) == a.data(k))
                        found = true;
                        break;
                      }
                  }

                  if (! found)
                  {
                    is_herm = false;
                    break;
                  }
                }
            }

            if (! diag_positive || ! is_herm)
            {
              is_herm = false;
              break;
            } 
          }

        if (is_herm)
          {
            if (typ == MatrixType::Full)
            typ = MatrixType::Hermitian;
            else if (typ == MatrixType::Banded)
            typ = MatrixType::Banded_Hermitian;
            else
            typ = MatrixType::Tridiagonal_Hermitian;
          }
      }
    }
}

MatrixType::MatrixType (const SparseComplexMatrix &a)
  : typ (MatrixType::Unknown),
    sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
    dense (false), full (false), nperm (0), perm (0)
{
  octave_idx_type nrows = a.rows ();
  octave_idx_type ncols = a.cols ();
  octave_idx_type nm = (ncols < nrows ? ncols : nrows);
  octave_idx_type nnz = a.nzmax ();

  if (octave_sparse_params::get_key ("spumoni") != 0.)
    (*current_liboctave_warning_handler) 
      ("Calculating Sparse Matrix Type");

  sp_bandden = octave_sparse_params::get_bandden();
  bool maybe_hermitian = false;
  typ = MatrixType::Full;

  if (nnz == nm)
    {
      matrix_type tmp_typ = MatrixType::Diagonal;
      octave_idx_type i;
      // Maybe the matrix is diagonal
      for (i = 0; i < nm; i++)
      {
        if (a.cidx(i+1) != a.cidx(i) + 1)
          {
            tmp_typ = MatrixType::Full;
            break;
          }
        if (a.ridx(i) != i)
          {
            tmp_typ = MatrixType::Permuted_Diagonal;
            break;
          }
      }
        
      if (tmp_typ == MatrixType::Permuted_Diagonal)
      {
        std::vector<bool> found (nrows);

        for (octave_idx_type j = 0; j < i; j++)
          found [j] = true;
        for (octave_idx_type j = i; j < nrows; j++)
          found [j] = false;
            
        for (octave_idx_type j = i; j < nm; j++)
          {
            if ((a.cidx(j+1) > a.cidx(j) + 1)  || 
              ((a.cidx(j+1) == a.cidx(j) + 1) && found [a.ridx(j)]))
            {
              tmp_typ = MatrixType::Full;
              break;
            }
            found [a.ridx(j)] = true;
          }
      }
      typ = tmp_typ;
    }

  if (typ == MatrixType::Full)
    {
      // Search for banded, upper and lower triangular matrices
      bool singular = false;
      upper_band = 0;
      lower_band = 0;
      for (octave_idx_type j = 0; j < ncols; j++)
      {
        bool zero_on_diagonal = false;
        if (j < nrows)
          {
            zero_on_diagonal = true;
            for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
            if (a.ridx(i) == j)
              {
                zero_on_diagonal = false;
                break;
              }
          }

        if (zero_on_diagonal)
          {
            singular = true;
            break;
          }

        if (a.cidx(j+1) != a.cidx(j))
          {
            octave_idx_type ru = a.ridx(a.cidx(j));
            octave_idx_type rl = a.ridx(a.cidx(j+1)-1);

            if (j - ru > upper_band)
            upper_band = j - ru;
              
            if (rl - j > lower_band)
            lower_band = rl - j;
          }
      }

      if (!singular)
      {
        bandden = double (nnz) /
          (double (ncols) * (double (lower_band) +
                         double (upper_band)) -
           0.5 * double (upper_band + 1) * double (upper_band) -
           0.5 * double (lower_band + 1) * double (lower_band));

        if (nrows == ncols && sp_bandden != 1. && bandden > sp_bandden)
          {
            if (upper_band == 1 && lower_band == 1)
            typ = MatrixType::Tridiagonal;
            else
            typ = MatrixType::Banded;

            octave_idx_type nnz_in_band = 
            (upper_band + lower_band + 1) * nrows -
            (1 + upper_band) * upper_band / 2 -
            (1 + lower_band) * lower_band / 2;
            if (nnz_in_band == nnz)
            dense = true;
            else 
            dense = false;
          }
        else if (upper_band == 0)
          typ = MatrixType::Lower;
        else if (lower_band == 0)
          typ = MatrixType::Upper;

        if (upper_band == lower_band && nrows == ncols)
          maybe_hermitian = true;
      }

      if (typ == MatrixType::Full)
      {
        // Search for a permuted triangular matrix, and test if
        // permutation is singular

        // FIXME
        // Perhaps this should be based on a dmperm algorithm
        bool found = false;

        nperm = ncols;
        perm = new octave_idx_type [ncols];

        for (octave_idx_type i = 0; i < ncols; i++)
          perm [i] = -1;

        for (octave_idx_type i = 0; i < nm; i++)
          {
            found = false;

            for (octave_idx_type j = 0; j < ncols; j++)
            {
              if ((a.cidx(j+1) - a.cidx(j)) > 0 && 
                  (a.ridx(a.cidx(j+1)-1) == i))
                {
                  perm [i] = j;
                  found = true;
                  break;
                }
            }

            if (!found)
            break;
          }

        if (found)
          {
            typ = MatrixType::Permuted_Upper;
            if (ncols > nrows)
            {
              octave_idx_type k = nrows;
              for (octave_idx_type i = 0; i < ncols; i++)
                if (perm [i] == -1)
                  perm[i] = k++;
            }
          }
        else if (a.cidx(nm) == a.cidx(ncols))
          {
            nperm = nrows;
            delete [] perm;
            perm = new octave_idx_type [nrows];
            OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, nrows);

            for (octave_idx_type i = 0; i < nrows; i++)
            {
              perm [i] = -1;
              tmp [i] = -1;
            }

            for (octave_idx_type j = 0; j < ncols; j++)
            for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
                perm [a.ridx(i)] = j;

            found = true;
            for (octave_idx_type i = 0; i < nm; i++)
            if (perm[i] == -1)
              {
                found = false;
                break;
              }
            else
              {
                tmp[perm[i]] = 1;
              }

            if (found)
            {
              octave_idx_type k = ncols;
              for (octave_idx_type i = 0; i < nrows; i++)
                {
                  if (tmp[i] == -1)
                  {
                    if (k < nrows)
                      {
                        perm[k++] = i;
                      }
                    else
                      {
                        found = false;
                        break;
                      }
                  }
                }
            }

            if (found)
            typ = MatrixType::Permuted_Lower;
            else
            {
              delete [] perm;
              nperm = 0;
            }
          }
        else
          {
            delete [] perm;
            nperm = 0;
          }
      }

      // FIXME
      // Disable lower under-determined and upper over-determined problems
      // as being detected, and force to treat as singular. As this seems
      // to cause issues
      if (((typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
         && nrows > ncols) ||
        ((typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
         && nrows < ncols))
      {
        typ = MatrixType::Rectangular;
        if (typ == MatrixType::Permuted_Upper ||
            typ == MatrixType::Permuted_Lower)
          delete [] perm;
        nperm = 0;
      }

      if (typ == MatrixType::Full && ncols != nrows)
      typ = MatrixType::Rectangular;

      if (maybe_hermitian && (typ == MatrixType::Full || 
                        typ == MatrixType::Tridiagonal || 
                        typ == MatrixType::Banded))
      {
        // Check for symmetry, with positive real diagonal, which
        // has a very good chance of being symmetric positive
        // definite..
        bool is_herm = true;

        for (octave_idx_type j = 0; j < ncols; j++)
          {
            bool diag_positive = false;

            for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
            {
              octave_idx_type ri = a.ridx(i);

              if (ri == j)
                {
                  if (a.data(i) == std::abs(a.data(i)))
                  diag_positive = true;
                  else
                  break;
                }
              else
                {
                  bool found = false;

                  for (octave_idx_type k = a.cidx(ri); k < a.cidx(ri+1); k++)
                  {
                    if (a.ridx(k) == j)
                      {
                        if (a.data(i) == conj(a.data(k)))
                        found = true;
                        break;
                      }
                  }

                  if (! found)
                  {
                    is_herm = false;
                    break;
                  }
                }
            }

            if (! diag_positive || ! is_herm)
            {
              is_herm = false;
              break;
            } 
          }

        if (is_herm)
          {
            if (typ == MatrixType::Full)
            typ = MatrixType::Hermitian;
            else if (typ == MatrixType::Banded)
            typ = MatrixType::Banded_Hermitian;
            else
            typ = MatrixType::Tridiagonal_Hermitian;
          }
      }
    }
}
MatrixType::MatrixType (const matrix_type t, bool _full)
  : typ (MatrixType::Unknown),
    sp_bandden (octave_sparse_params::get_bandden()),
    bandden (0), upper_band (0), lower_band (0),
    dense (false), full (_full), nperm (0), perm (0)
{
  if (t == MatrixType::Full || t == MatrixType::Diagonal ||
      t == MatrixType::Permuted_Diagonal || t == MatrixType::Upper ||
      t == MatrixType::Lower || t == MatrixType::Tridiagonal ||
      t == MatrixType::Tridiagonal_Hermitian || t == MatrixType::Rectangular)
    typ = t;
  else
    (*current_liboctave_warning_handler) ("Invalid matrix type");
}

MatrixType::MatrixType (const matrix_type t, const octave_idx_type np,
                  const octave_idx_type *p, bool _full)
  : typ (MatrixType::Unknown),
    sp_bandden (octave_sparse_params::get_bandden()),
    bandden (0), upper_band (0), lower_band (0),
    dense (false), full (_full), nperm (0), perm (0)
{
  if ((t == MatrixType::Permuted_Upper || t == MatrixType::Permuted_Lower) &&
      np > 0 && p != 0)
    {
      typ = t;
      nperm = np;
      perm = new octave_idx_type [nperm];
      for (octave_idx_type i = 0; i < nperm; i++)
      perm[i] = p[i];
    }
  else
    (*current_liboctave_warning_handler) ("Invalid matrix type");
}

MatrixType::MatrixType (const matrix_type t, const octave_idx_type ku,
                  const octave_idx_type kl, bool _full)
  : typ (MatrixType::Unknown),
    sp_bandden (octave_sparse_params::get_bandden()),
    bandden (0), upper_band (0), lower_band (0),
    dense (false), full (_full), nperm (0), perm (0)
{
  if (t == MatrixType::Banded || t == MatrixType::Banded_Hermitian)
    {
      typ = t;
      upper_band = ku;
      lower_band = kl;
    }
  else
    (*current_liboctave_warning_handler) ("Invalid sparse matrix type"); 
}

MatrixType::~MatrixType (void) 
{ 
  if (nperm != 0)
    {
      delete [] perm; 
    }
}

MatrixType& 
MatrixType::operator = (const MatrixType& a)
{
  if (this != &a)
    {
      typ = a.typ;
      sp_bandden = a.sp_bandden;
      bandden = a.bandden;
      upper_band = a.upper_band;
      lower_band = a.lower_band;
      dense = a.dense;
      full = a.full;
      nperm = a.nperm;

      if (nperm != 0)
      {
        perm = new octave_idx_type [nperm];
        for (octave_idx_type i = 0; i < nperm; i++)
          perm[i] = a.perm[i];
      }
    }

  return *this;
}

int
MatrixType::type (bool quiet)
{
  if (typ != MatrixType::Unknown && (full ||
      sp_bandden == octave_sparse_params::get_bandden()))
    {
      if (!quiet &&
        octave_sparse_params::get_key ("spumoni") != 0.)
      (*current_liboctave_warning_handler) 
        ("Using Cached Matrix Type");
      
      return typ;
    }

  if (typ != MatrixType::Unknown && 
      octave_sparse_params::get_key ("spumoni") != 0.)
    (*current_liboctave_warning_handler) 
      ("Invalidating Matrix Type");

  typ = MatrixType::Unknown;

  return typ;
}

int
MatrixType::type (const SparseMatrix &a)
{
  if (typ != MatrixType::Unknown && (full ||
      sp_bandden == octave_sparse_params::get_bandden()))
    {
      if (octave_sparse_params::get_key ("spumoni") != 0.)
      (*current_liboctave_warning_handler) 
        ("Using Cached Matrix Type");
      
      return typ;
    }

  MatrixType tmp_typ (a);
  typ = tmp_typ.typ;
  sp_bandden = tmp_typ.sp_bandden;
  bandden = tmp_typ.bandden;
  upper_band = tmp_typ.upper_band;
  lower_band = tmp_typ.lower_band;
  dense = tmp_typ.dense;
  full = tmp_typ.full;
  nperm = tmp_typ.nperm;

  if (nperm != 0)
    {
      perm = new octave_idx_type [nperm];
      for (octave_idx_type i = 0; i < nperm; i++)
      perm[i] = tmp_typ.perm[i];
    }

  return typ;
}

int
MatrixType::type (const SparseComplexMatrix &a)
{
  if (typ != MatrixType::Unknown && (full || 
      sp_bandden == octave_sparse_params::get_bandden()))
    {
      if (octave_sparse_params::get_key ("spumoni") != 0.)
      (*current_liboctave_warning_handler) 
        ("Using Cached Matrix Type");
      
      return typ;
    }

  MatrixType tmp_typ (a);
  typ = tmp_typ.typ;
  sp_bandden = tmp_typ.sp_bandden;
  bandden = tmp_typ.bandden;
  upper_band = tmp_typ.upper_band;
  lower_band = tmp_typ.lower_band;
  dense = tmp_typ.dense;
  full = tmp_typ.full;
  nperm = tmp_typ.nperm;

  if (nperm != 0)
    {
      perm = new octave_idx_type [nperm];
      for (octave_idx_type i = 0; i < nperm; i++)
      perm[i] = tmp_typ.perm[i];
    }

  return typ;
}
int
MatrixType::type (const Matrix &a)
{
  if (typ != MatrixType::Unknown)
    {
      if (octave_sparse_params::get_key ("spumoni") != 0.)
      (*current_liboctave_warning_handler) 
        ("Using Cached Matrix Type");
      
      return typ;
    }

  MatrixType tmp_typ (a);
  typ = tmp_typ.typ;
  full = tmp_typ.full;
  nperm = tmp_typ.nperm;

  if (nperm != 0)
    {
      perm = new octave_idx_type [nperm];
      for (octave_idx_type i = 0; i < nperm; i++)
      perm[i] = tmp_typ.perm[i];
    }

  return typ;
}

int
MatrixType::type (const ComplexMatrix &a)
{
  if (typ != MatrixType::Unknown)
    {
      if (octave_sparse_params::get_key ("spumoni") != 0.)
      (*current_liboctave_warning_handler) 
        ("Using Cached Matrix Type");
      
      return typ;
    }

  MatrixType tmp_typ (a);
  typ = tmp_typ.typ;
  full = tmp_typ.full; 
  nperm = tmp_typ.nperm;

  if (nperm != 0)
    {
      perm = new octave_idx_type [nperm];
      for (octave_idx_type i = 0; i < nperm; i++)
      perm[i] = tmp_typ.perm[i];
    }

  return typ;
}

void
MatrixType::info () const
{
  if (octave_sparse_params::get_key ("spumoni") != 0.)
    {
      if (typ == MatrixType::Unknown)
      (*current_liboctave_warning_handler) 
        ("Unknown Matrix Type");
      else if (typ == MatrixType::Diagonal)
      (*current_liboctave_warning_handler) 
        ("Diagonal Sparse Matrix");
      else if (typ == MatrixType::Permuted_Diagonal)
      (*current_liboctave_warning_handler) 
        ("Permuted Diagonal Sparse Matrix");
      else if (typ == MatrixType::Upper)
      (*current_liboctave_warning_handler) 
        ("Upper Triangular Matrix");
      else if (typ == MatrixType::Lower)
      (*current_liboctave_warning_handler) 
        ("Lower Triangular Matrix");
      else if (typ == MatrixType::Permuted_Upper)
      (*current_liboctave_warning_handler) 
        ("Permuted Upper Triangular Matrix");
      else if (typ == MatrixType::Permuted_Lower)
      (*current_liboctave_warning_handler) 
        ("Permuted Lower Triangular Matrix");
      else if (typ == MatrixType::Banded)
      (*current_liboctave_warning_handler) 
        ("Banded Sparse Matrix %d-1-%d (Density %f)", lower_band, 
         upper_band, bandden);
      else if (typ == MatrixType::Banded_Hermitian)
      (*current_liboctave_warning_handler) 
        ("Banded Hermitian/Symmetric Sparse Matrix %d-1-%d (Density %f)", 
         lower_band, upper_band, bandden);
      else if (typ == MatrixType::Hermitian)
      (*current_liboctave_warning_handler) 
        ("Hermitian/Symmetric Matrix");
      else if (typ == MatrixType::Tridiagonal)
      (*current_liboctave_warning_handler) 
        ("Tridiagonal Sparse Matrix");
      else if (typ == MatrixType::Tridiagonal_Hermitian)
      (*current_liboctave_warning_handler) 
        ("Hermitian/Symmetric Tridiagonal Sparse Matrix");
      else if (typ == MatrixType::Rectangular)
      (*current_liboctave_warning_handler) 
        ("Rectangular/Singular Matrix");
      else if (typ == MatrixType::Full)
      (*current_liboctave_warning_handler) 
        ("Full Matrix");
    }
}

void
MatrixType::mark_as_symmetric (void)
{
  if (typ == MatrixType::Tridiagonal || 
      typ == MatrixType::Tridiagonal_Hermitian)
    typ = MatrixType::Tridiagonal_Hermitian;
  else if (typ == MatrixType::Banded ||
         typ == MatrixType::Banded_Hermitian)
    typ = MatrixType::Banded_Hermitian;
  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian || 
         typ == MatrixType::Unknown)
    typ = MatrixType::Hermitian;
  else
    (*current_liboctave_error_handler) 
      ("Can not mark current matrix type as symmetric");
}

void
MatrixType::mark_as_unsymmetric (void)
{
  if (typ == MatrixType::Tridiagonal || 
      typ == MatrixType::Tridiagonal_Hermitian)
    typ = MatrixType::Tridiagonal;
  else if (typ == MatrixType::Banded ||
         typ == MatrixType::Banded_Hermitian)
    typ = MatrixType::Banded;
  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian || 
         typ == MatrixType::Unknown)
    typ = MatrixType::Full;
}

void
MatrixType::mark_as_permuted (const octave_idx_type np, const octave_idx_type *p)
{
  nperm = np;
  perm = new octave_idx_type [nperm];
  for (octave_idx_type i = 0; i < nperm; i++)
    perm[i] = p[i];

  if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
    typ = MatrixType::Permuted_Diagonal;
  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
    typ = MatrixType::Permuted_Upper;
  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
    typ = MatrixType::Permuted_Lower;
  else
    (*current_liboctave_error_handler) 
      ("Can not mark current matrix type as symmetric");
}

void
MatrixType::mark_as_unpermuted (void)
{
  if (nperm)
    {
      nperm = 0;
      delete [] perm;
    }

  if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
    typ = MatrixType::Diagonal;
  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
    typ = MatrixType::Upper;
  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
    typ = MatrixType::Lower;
}

MatrixType
MatrixType::transpose (void) const
{
  MatrixType retval (*this);
  if (typ == MatrixType::Upper)
    retval.typ = MatrixType::Lower;
  else if (typ == MatrixType::Permuted_Upper)
    retval.typ = MatrixType::Permuted_Lower;
  else if (typ == MatrixType::Lower)
    retval.typ = MatrixType::Upper;
  else if (typ == MatrixType::Permuted_Lower)
    retval.typ = MatrixType::Permuted_Upper;
  else if (typ == MatrixType::Banded)
    {
      retval.upper_band = lower_band;
      retval.lower_band = upper_band;
    }

  return retval;
}

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/


Generated by  Doxygen 1.6.0   Back to index