matrix.cpp

Go to the documentation of this file.

/*
 * matrix.cpp - matrix class implementation
 *
 * Copyright (C) 2003-2009 Stefan Jahn <stefan@lkcc.org>
 *
 * This 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.
 *
 * This software 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 package; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
 * Boston, MA 02110-1301, USA.
 *
 * $Id$
 *
 */
#if HAVE_CONFIG_H
# include <config.h>
#endif

#include <assert.h>
#include <stdio.h>
#include <cstdlib>
#include <string.h>
#include <cmath>

#include "logging.h"
#include "object.h"
#include "complex.h"
#include "vector.h"
#include "matrix.h"

namespace qucs {

matrix::matrix () {
  rows = 0;
  cols = 0;
  data = NULL;
}

matrix::matrix (int s)  {
  rows = cols = s;
  data = (s > 0) ? new nr_complex_t[s * s] : NULL;
}

/* \brief Creates a matrix

   Constructor creates an unnamed instance of the matrix class with a
   certain number of rows and columns.
   \param[in] r number of rows
   \param[in] c number of column
   \todo Why not r and c constant
   \todo Assert r >= 0 and c >= 0
*/
matrix::matrix (int r, int c)  {
  rows = r;
  cols = c;
  data = (r > 0 && c > 0) ? new nr_complex_t[r * c] : NULL;
}

/* \brief copy constructor

   The copy constructor creates a new instance based on the given
   matrix object.
   \todo Add assert tests
*/
matrix::matrix (const matrix & m) {
  rows = m.rows;
  cols = m.cols;
  data = NULL;

  // copy matrix elements
  if (rows > 0 && cols > 0) {
    data = new nr_complex_t[rows * cols];
    memcpy (data, m.data, sizeof (nr_complex_t) * rows * cols);
  }
}

const matrix& matrix::operator=(const matrix & m) {
  if (&m != this) {
    rows = m.rows;
    cols = m.cols;
    if (data) {
      delete[] data;
      data = NULL;
    }
    if (rows > 0 && cols > 0) {
      data = new nr_complex_t[rows * cols];
      memcpy (data, m.data, sizeof (nr_complex_t) * rows * cols);
    }
  }
  return *this;
}

matrix::~matrix () {
  if (data) delete[] data;
}

nr_complex_t matrix::get (int r, int c) {
  return data[r * cols + c];
}

void matrix::set (int r, int c, nr_complex_t z) {
  data[r * cols + c] = z;
}

#ifdef DEBUG

void matrix::print (void) {
  for (int r = 0; r < rows; r++) {
    for (int c = 0; c < cols; c++) {
      fprintf (stderr, "%+.2e,%+.2e ", (double) real (get (r, c)),
               (double) imag (get (r, c)));
    }
    fprintf (stderr, "\n");
  }
}
#endif /* DEBUG */

matrix operator + (matrix a, matrix b) {
  assert (a.getRows () == b.getRows () && a.getCols () == b.getCols ());

  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, a.get (r, c) + b.get (r, c));
  return res;
}

matrix matrix::operator += (matrix a) {
  assert (a.getRows () == rows && a.getCols () == cols);

  int r, c, i;
  for (i = 0, r = 0; r < a.getRows (); r++)
    for (c = 0; c < a.getCols (); c++, i++)
      data[i] += a.get (r, c);
  return *this;
}

matrix operator - (matrix a, matrix b) {
  assert (a.getRows () == b.getRows () && a.getCols () == b.getCols ());

  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, a.get (r, c) - b.get (r, c));
  return res;
}

matrix matrix::operator - () {
  matrix res (getRows (), getCols ());
  int r, c, i;
  for (i = 0, r = 0; r < getRows (); r++)
    for (c = 0; c < getCols (); c++, i++)
      res.set (r, c, -data[i]);
  return res;
}

matrix matrix::operator -= (matrix a) {
  assert (a.getRows () == rows && a.getCols () == cols);
  int r, c, i;
  for (i = 0, r = 0; r < a.getRows (); r++)
    for (c = 0; c < a.getCols (); c++, i++)
      data[i] -= a.get (r, c);
  return *this;
}

matrix operator * (matrix a, nr_complex_t z) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, a.get (r, c) * z);
  return res;
}

matrix operator * (nr_complex_t z, matrix a) {
  return a * z;
}

matrix operator * (matrix a, nr_double_t d) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, a.get (r, c) * d);
  return res;
}

matrix operator * (nr_double_t d, matrix a) {
  return a * d;
}

matrix operator / (matrix a, nr_complex_t z) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, a.get (r, c) / z);
  return res;
}

matrix operator / (matrix a, nr_double_t d) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, a.get (r, c) / d);
  return res;
}

matrix operator * (matrix a, matrix b) {
  assert (a.getCols () == b.getRows ());

  int r, c, i, n = a.getCols ();
  nr_complex_t z;
  matrix res (a.getRows (), b.getCols ());
  for (r = 0; r < a.getRows (); r++) {
    for (c = 0; c < b.getCols (); c++) {
      for (i = 0, z = 0; i < n; i++) z += a.get (r, i) * b.get (i, c);
      res.set (r, c, z);
    }
  }
  return res;
}

matrix operator + (matrix a, nr_complex_t z) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, a.get (r, c) + z);
  return res;
}

matrix operator + (nr_complex_t z, matrix a) {
  return a + z;
}

matrix operator + (matrix a, nr_double_t d) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, a.get (r, c) + d);
  return res;
}

matrix operator + (nr_double_t d, matrix a) {
  return a + d;
}

matrix operator - (matrix a, nr_complex_t z) {
  return -z + a;
}

matrix operator - (nr_complex_t z, matrix a) {
  return -a + z;
}

matrix operator - (matrix a, nr_double_t d) {
  return -d + a;
}

matrix operator - (nr_double_t d, matrix a) {
  return -a + d;
}

matrix transpose (matrix a) {
  matrix res (a.getCols (), a.getRows ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (c, r, a.get (r, c));
  return res;
}

matrix conj (matrix a) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, conj (a.get (r, c)));
  return res;
}

matrix adjoint (matrix a) {
  return transpose (conj (a));
}

matrix abs (matrix a) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, abs (a.get (r, c)));
  return res;
}

matrix dB (matrix a) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, dB (a.get (r, c)));
  return res;
}

matrix arg (matrix a) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, arg (a.get (r, c)));
  return res;
}

matrix real (matrix a) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, real (a.get (r, c)));
  return res;
}

matrix imag (matrix a) {
  matrix res (a.getRows (), a.getCols ());
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, imag (a.get (r, c)));
  return res;
}

matrix sqr (matrix a) {
  return a * a;
}

matrix eye (int rs, int cs) {
  matrix res (rs, cs);
  for (int r = 0; r < res.getRows (); r++)
    for (int c = 0; c < res.getCols (); c++)
      if (r == c) res.set (r, c, 1);
  return res;
}

matrix eye (int s) {
  return eye (s, s);
}

matrix diagonal (qucs::vector diag) {
  int size = diag.getSize ();
  matrix res (size);
  for (int i = 0; i < size; i++) res (i, i) = diag (i);
  return res;
}

// Compute n-th power of the given matrix.
matrix pow (matrix a, int n) {
  matrix res;
  if (n == 0) {
    res = eye (a.getRows (), a.getCols ());
  }
  else {
    res = a = n < 0 ? inverse (a) : a;
    for (int i = 1; i < std::abs (n); i++)
      res = res * a;
  }
  return res;
}

nr_complex_t cofactor (matrix a, int u, int v) {
  matrix res (a.getRows () - 1, a.getCols () - 1);
  int r, c, ra, ca;
  for (ra = r = 0; r < res.getRows (); r++, ra++) {
    if (ra == u) ra++;
    for (ca = c = 0; c < res.getCols (); c++, ca++) {
      if (ca == v) ca++;
      res.set (r, c, a.get (ra, ca));
    }
  }
  nr_complex_t z = detLaplace (res);
  return ((u + v) & 1) ? -z : z;
}

nr_complex_t detLaplace (matrix a) {
  assert (a.getRows () == a.getCols ());
  int s = a.getRows ();
  nr_complex_t res = 0;
  if (s > 1) {
    /* always use the first row for sub-determinant, but you should
       use the row or column with most zeros in it */
    int r = 0;
    for (int i = 0; i < s; i++) {
      res += a.get (r, i) * cofactor (a, r, i);
    }
    return res;
  }
  /* 1 by 1 matrix */
  else if (s == 1) {
    return a (0, 0);
  }
  /* 0 by 0 matrix */
  return 1;
}

nr_complex_t detGauss (matrix a) {
  assert (a.getRows () == a.getCols ());
  nr_double_t MaxPivot;
  nr_complex_t f, res;
  matrix b;
  int i, c, r, pivot, n = a.getCols ();

  // return special matrix cases
  if (n == 0) return 1;
  if (n == 1) return a (0, 0);

  // make copy of original matrix
  b = matrix (a);

  // triangulate the matrix
  for (res = 1, i = 0; i < n; i++) {
    // find maximum column value for pivoting
    for (MaxPivot = 0, pivot = r = i; r < n; r++) {
      if (abs (b.get (r, i)) > MaxPivot) {
        MaxPivot = abs (b.get (r, i));
        pivot = r;
      }
    }
    // exchange rows if necessary
    assert (MaxPivot != 0);
    if (i != pivot) { b.exchangeRows (i, pivot); res = -res; }
    // compute new rows and columns
    for (r = i + 1; r < n; r++) {
      f = b.get (r, i) / b.get (i, i);
      for (c = i + 1; c < n; c++) {
        b.set (r, c, b.get (r, c) - f * b.get (i, c));
      }
    }
  }

  // now compute determinant by multiplying diagonal
  for (i = 0; i < n; i++) res *= b.get (i, i);
  return res;
}

nr_complex_t det (matrix a) {
#if 0
  return detLaplace (a);
#else
  return detGauss (a);
#endif
}

matrix inverseLaplace (matrix a) {
  matrix res (a.getRows (), a.getCols ());
  nr_complex_t d = detLaplace (a);
  assert (abs (d) != 0); // singular matrix
  for (int r = 0; r < a.getRows (); r++)
    for (int c = 0; c < a.getCols (); c++)
      res.set (r, c, cofactor (a, c, r) / d);
  return res;
}

matrix inverseGaussJordan (matrix a) {
  nr_double_t MaxPivot;
  nr_complex_t f;
  matrix b, e;
  int i, c, r, pivot, n = a.getCols ();

  // create temporary matrix and the result matrix
  b = matrix (a);
  e = eye (n);

  // create the eye matrix in 'b' and the result in 'e'
  for (i = 0; i < n; i++) {
    // find maximum column value for pivoting
    for (MaxPivot = 0, pivot = r = i; r < n; r++) {
      if (abs (b (r, i)) > MaxPivot) {
        MaxPivot = abs (b (r, i));
        pivot = r;
      }
    }
    // exchange rows if necessary
    assert (MaxPivot != 0); // singular matrix
    if (i != pivot) {
      b.exchangeRows (i, pivot);
      e.exchangeRows (i, pivot);
    }

    // compute current row
    for (f = b (i, i), c = 0; c < n; c++) {
      b (i, c) /= f;
      e (i, c) /= f;
    }

    // compute new rows and columns
    for (r = 0; r < n; r++) {
      if (r != i) {
        for (f = b (r, i), c = 0; c < n; c++) {
          b (r, c) -= f * b (i, c);
          e (r, c) -= f * e (i, c);
        }
      }
    }
  }
  return e;
}

matrix inverse (matrix a) {
#if 0
  return inverseLaplace (a);
#else
  return inverseGaussJordan (a);
#endif
}

matrix stos (matrix s, qucs::vector zref, qucs::vector z0) {
  int d = s.getRows ();
  matrix e, r, a;

  assert (d == s.getCols () && d == z0.getSize () && d == zref.getSize ());

  e = eye (d);
  r = diagonal ((z0 - zref) / (z0 + zref));
  a = diagonal (sqrt (z0 / zref) / (z0 + zref));
  return inverse (a) * (s - r) * inverse (e - r * s) * a;
}

matrix stos (matrix s, nr_complex_t zref, nr_complex_t z0) {
  int d = s.getRows ();
  return stos (s, qucs::vector (d, zref), qucs::vector (d, z0));
}

matrix stos (matrix s, nr_double_t zref, nr_double_t z0) {
  return stos (s, nr_complex_t (zref, 0), nr_complex_t (z0, 0));
}

matrix stos (matrix s, qucs::vector zref, nr_complex_t z0) {
  return stos (s, zref, qucs::vector (zref.getSize (), z0));
}

matrix stos (matrix s, nr_complex_t zref, qucs::vector z0) {
  return stos (s, qucs::vector (z0.getSize (), zref), z0);
}

matrix stoz (matrix s, qucs::vector z0) {
  int d = s.getRows ();
  matrix e, zref, gref;

  assert (d == s.getCols () && d == z0.getSize ());

  e = eye (d);
  zref = diagonal (z0);
  gref = diagonal (sqrt (real (1 / z0)));
  return inverse (gref) * inverse (e - s) * (s * zref + zref) * gref;
}

matrix stoz (matrix s, nr_complex_t z0) {
  return stoz (s, qucs::vector (s.getRows (), z0));
}

matrix ztos (matrix z, qucs::vector z0) {
  int d = z.getRows ();
  matrix e, zref, gref;

  assert (d == z.getCols () && d == z0.getSize ());

  e = eye (d);
  zref = diagonal (z0);
  gref = diagonal (sqrt (real (1 / z0)));
  return gref * (z - zref) * inverse (z + zref) * inverse (gref);
}

matrix ztos (matrix z, nr_complex_t z0) {
  return ztos (z, qucs::vector (z.getRows (), z0));
}

matrix ztoy (matrix z) {
  assert (z.getRows () == z.getCols ());
  return inverse (z);
}

matrix stoy (matrix s, qucs::vector z0) {
  int d = s.getRows ();
  matrix e, zref, gref;

  assert (d == s.getCols () && d == z0.getSize ());

  e = eye (d);
  zref = diagonal (z0);
  gref = diagonal (sqrt (real (1 / z0)));
  return inverse (gref) * inverse (s * zref + zref) * (e - s) * gref;
}

matrix stoy (matrix s, nr_complex_t z0) {
  return stoy (s, qucs::vector (s.getRows (), z0));
}

matrix ytos (matrix y, qucs::vector z0) {
  int d = y.getRows ();
  matrix e, zref, gref;

  assert (d == y.getCols () && d == z0.getSize ());

  e = eye (d);
  zref = diagonal (z0);
  gref = diagonal (sqrt (real (1 / z0)));
  return gref * (e - zref * y) * inverse (e + zref * y) * inverse (gref);
}
matrix ytos (matrix y, nr_complex_t z0) {
  return ytos (y, qucs::vector (y.getRows (), z0));
}
matrix stoa (matrix s, nr_complex_t z1, nr_complex_t z2) {
  nr_complex_t d = s (0, 0) * s (1, 1) - s (0, 1) * s (1, 0);
  nr_complex_t n = 2.0 * s (1, 0) * sqrt (fabs (real (z1) * real (z2)));
  matrix a (2);

  assert (s.getRows () >= 2 && s.getCols () >= 2);

  a.set (0, 0, (conj (z1) + z1 * s (0, 0) -
                conj (z1) * s (1, 1) - z1 * d) / n);
  a.set (0, 1, (conj (z1) * conj (z2) + z1 * conj (z2) * s (0, 0) +
                conj (z1) * z2 * s (1, 1) + z1 * z2 * d) / n);
  a.set (1, 0, (1.0 - s (0, 0) - s (1, 1) + d) / n);
  a.set (1, 1, (conj (z2) - conj (z2) * s (0, 0) +
                z2 * s (1, 1) - z2 * d) / n);
  return a;
}


matrix atos (matrix a, nr_complex_t z1, nr_complex_t z2) {
  nr_complex_t d = 2.0 * sqrt (fabs (real (z1) * real (z2)));
  nr_complex_t n = a (0, 0) * z2 + a (0, 1) +
    a (1, 0) * z1 * z2 + a (1, 1) * z1;
  matrix s (2);

  assert (a.getRows () >= 2 && a.getCols () >= 2);

  s.set (0, 0, (a (0, 0) * z2 + a (0, 1)
                - a (1, 0) * conj (z1) * z2 - a (1, 1) * conj (z1)) / n);
  s.set (0, 1, (a (0, 0) * a (1, 1) -
                a (0, 1) * a (1, 0)) * d / n);
  s.set (1, 0, d / n);
  s.set (1, 1, (a (1, 1) * z1 - a (0, 0) * conj (z2) +
                a (0, 1) - a (1, 0) * z1 * conj (z2)) / n);
  return s;
}

matrix stoh (matrix s, nr_complex_t z1, nr_complex_t z2) {
  nr_complex_t n = s (0, 1) * s (1, 0);
  nr_complex_t d = (1.0 - s (0, 0)) * (1.0 + s (1, 1)) + n;
  matrix h (2);

  assert (s.getRows () >= 2 && s.getCols () >= 2);

  h.set (0, 0, ((1.0 + s (0, 0)) * (1.0 + s (1, 1)) - n) * z1 / d);
  h.set (0, 1, +2.0 * s (0, 1) / d);
  h.set (1, 0, -2.0 * s (1, 0) / d);
  h.set (1, 1, ((1.0 - s (0, 0)) * (1.0 - s (1, 1)) - n) / z2 / d);
  return h;
}

matrix htos (matrix h, nr_complex_t z1, nr_complex_t z2) {
  nr_complex_t n = h (0, 1) * h (1, 0);
  nr_complex_t d = (1.0 + h (0, 0) / z1) * (1.0 + z2 * h (1, 1)) - n;
  matrix s (2);

  assert (h.getRows () >= 2 && h.getCols () >= 2);

  s.set (0, 0, ((h (0, 0) / z1 - 1.0) * (1.0 + z2 * h (1, 1)) - n) / d);
  s.set (0, 1, +2.0 * h (0, 1) / d);
  s.set (1, 0, -2.0 * h (1, 0) / d);
  s.set (1, 1, ((1.0 + h (0, 0) / z1) * (1.0 - z2 * h (1, 1)) + n) / d);
  return s;
}

/*\brief Converts scattering parameters to second hybrid matrix.
  \bug{Programmed formulae are valid only for Z real}
  \bug{Not documented and references}
  \param[in] s Scattering matrix
  \param[in] z1 impedance at input 1
  \param[in] z2 impedance at input 2
  \return second hybrid matrix
  \note Assert 2 by 2 matrix
  \todo Why not s,z1,z2 const
*/
matrix stog (matrix s, nr_complex_t z1, nr_complex_t z2) {
  nr_complex_t n = s (0, 1) * s (1, 0);
  nr_complex_t d = (1.0 + s (0, 0)) * (1.0 - s (1, 1)) + n;
  matrix g (2);

  assert (s.getRows () >= 2 && s.getCols () >= 2);

  g.set (0, 0, ((1.0 - s (0, 0)) * (1.0 - s (1, 1)) - n) / z1 / d);
  g.set (0, 1, -2.0 * s (0, 1) / d);
  g.set (1, 0, +2.0 * s (1, 0) / d);
  g.set (1, 1, ((1.0 + s (0, 0)) * (1.0 + s (1, 1)) - n) * z2 / d);
  return g;
}

/*\brief Converts second hybrid matrix to scattering parameters.
  \bug{Programmed formulae are valid only for Z real}
  \bug{Not documented and references}
  \param[in] g second hybrid matrix
  \param[in] z1 impedance at input 1
  \param[in] z2 impedance at input 2
  \return scattering matrix
  \note Assert 2 by 2 matrix
  \todo Why not g,z1,z2 const
*/
matrix gtos (matrix g, nr_complex_t z1, nr_complex_t z2) {
  nr_complex_t n = g (0, 1) * g (1, 0);
  nr_complex_t d = (1.0 + g (0, 0) * z1) * (1.0 + g (1, 1) / z2) - n;
  matrix s (2);

  assert (g.getRows () >= 2 && g.getCols () >= 2);

  s.set (0, 0, ((1.0 - g (0, 0) * z1) * (1.0 + g (1, 1) / z2) + n) / d);
  s.set (0, 1, -2.0 * g (0, 1) / d);
  s.set (1, 0, +2.0 * g (1, 0) / d);
  s.set (1, 1, ((g (0, 0) * z1 + 1.0) * (g (1, 1) / z2 - 1.0) - n) / d);
  return s;
}

matrix ytoz (matrix y) {
  assert (y.getRows () == y.getCols ());
  return inverse (y);
}

matrix cytocs (matrix cy, matrix s) {
  matrix e = eye (s.getRows ());

  assert (cy.getRows () == cy.getCols () && s.getRows () == s.getCols () &&
          cy.getRows () == s.getRows ());

  return (e + s) * cy * adjoint (e + s) / 4;
}

matrix cstocy (matrix cs, matrix y) {
  matrix e = eye (y.getRows ());

  assert (cs.getRows () == cs.getCols () && y.getRows () == y.getCols () &&
          cs.getRows () == y.getRows ());

  return (e + y) * cs * adjoint (e + y);
}

matrix cztocs (matrix cz, matrix s) {
  matrix e = eye (s.getRows ());

  assert (cz.getRows () == cz.getCols () && s.getRows () == s.getCols () &&
          cz.getRows () == s.getRows ());

  return (e - s) * cz * adjoint (e - s) / 4;
}

matrix cstocz (matrix cs, matrix z) {
  assert (cs.getRows () == cs.getCols () && z.getRows () == z.getCols () &&
          cs.getRows () == z.getRows ());
  matrix e = eye (z.getRows ());
  return (e + z) * cs * adjoint (e + z);
}

matrix cztocy (matrix cz, matrix y) {
  assert (cz.getRows () == cz.getCols () && y.getRows () == y.getCols () &&
          cz.getRows () == y.getRows ());

  return y * cz * adjoint (y);
}

matrix cytocz (matrix cy, matrix z) {
  assert (cy.getRows () == cy.getCols () && z.getRows () == z.getCols () &&
          cy.getRows () == z.getRows ());
  return z * cy * adjoint (z);
}

void matrix::exchangeRows (int r1, int r2) {
  nr_complex_t * s = new nr_complex_t[cols];
  int len = sizeof (nr_complex_t) * cols;

  assert (r1 >= 0 && r2 >= 0 && r1 < rows && r2 < rows);

  memcpy (s, &data[r1 * cols], len);
  memcpy (&data[r1 * cols], &data[r2 * cols], len);
  memcpy (&data[r2 * cols], s, len);
  delete[] s;
}

void matrix::exchangeCols (int c1, int c2) {
  nr_complex_t s;

  assert (c1 >= 0 && c2 >= 0 && c1 < cols && c2 < cols);

  for (int r = 0; r < rows * cols; r += cols) {
    s = data[r + c1];
    data[r + c1] = data[r + c2];
    data[r + c2] = s;
  }
}

matrix twoport (matrix m, char in, char out) {
  assert (m.getRows () >= 2 && m.getCols () >= 2);
  nr_complex_t d;
  matrix res (2);

  switch (in) {
  case 'Y':
    switch (out) {
    case 'Y': // Y to Y
      res = m;
      break;
    case 'Z': // Y to Z
      d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
      res.set (0, 0, m (1, 1) / d);
      res.set (0, 1, -m (0, 1) / d);
      res.set (1, 0, -m (1, 0) / d);
      res.set (1, 1, m (0, 0) / d);
      break;
    case 'H': // Y to H
      d = m (0, 0);
      res.set (0, 0, 1.0 / d);
      res.set (0, 1, -m (0, 1) / d);
      res.set (1, 0, m (1, 0) / d);
      res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d);
      break;
    case 'G': // Y to G
      d = m (1, 1);
      res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
      res.set (0, 1, m (0, 1) / d);
      res.set (1, 0, -m (1, 0) / d);
      res.set (1, 1, 1.0 / d);
      break;
    case 'A': // Y to A
      d = m (1, 0);
      res.set (0, 0, -m (1, 1) / d);
      res.set (0, 1, -1.0 / d);
      res.set (1, 0, m (0, 1) - m (1, 1) * m (0, 0) / d);
      res.set (1, 1, -m (0, 0) / d);
      break;
    case 'S': // Y to S
      res = ytos (m);
      break;
    }
    break;
  case 'Z':
    switch (out) {
    case 'Y': // Z to Y
      d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
      res.set (0, 0, m (1, 1) / d);
      res.set (0, 1, -m (0, 1) / d);
      res.set (1, 0, -m (1, 0) / d);
      res.set (1, 1, m (0, 0) / d);
      break;
    case 'Z': // Z to Z
      res = m;
      break;
    case 'H': // Z to H
      d = m (1, 1);
      res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
      res.set (0, 1, m (0, 1) / d);
      res.set (1, 0, -m (1, 0) / d);
      res.set (1, 1, 1.0 / d);
      break;
    case 'G': // Z to G
      d = m (0, 0);
      res.set (0, 0, 1.0 / d);
      res.set (0, 1, -m (0, 1) / d);
      res.set (1, 0, m (1, 0) / d);
      res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d);
      break;
    case 'A': // Z to A
      d = m (1, 0);
      res.set (0, 0, m (0, 0) / d);
      res.set (0, 1, m (0, 0) * m (1, 1) / d - m (0, 1));
      res.set (1, 0, 1.0 / d);
      res.set (1, 1, m (1, 1) / d);
      break;
    case 'S': // Z to S
      res = ztos (m);
      break;
    }
    break;
  case 'H':
    switch (out) {
    case 'Y': // H to Y
      d = m (0, 0);
      res.set (0, 0, 1.0 / d);
      res.set (0, 1, -m (0, 1) / d);
      res.set (1, 0, m (1, 0) / d);
      res.set (1, 1, m (1, 1) - m (0, 1) * m.get(2, 1) / d);
      break;
    case 'Z': // H to Z
      d = m (1, 1);
      res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
      res.set (0, 1, m (0, 1) / d);
      res.set (1, 0, -m (1, 0) / d);
      res.set (1, 1, 1.0 / d);
      break;
    case 'H': // H to H
      res = m;
      break;
    case 'G': // H to G
      d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
      res.set (0, 0, m (1, 1) / d);
      res.set (0, 1, -m (0, 1) / d);
      res.set (1, 0, -m (1, 0) / d);
      res.set (1, 1, m (0, 0) / d);
      break;
    case 'A': // H to A
      d = m (1, 0);
      res.set (0, 0, m (0, 1) - m (0, 0) * m (1, 1) / d);
      res.set (0, 1, -m (0, 0) / d);
      res.set (1, 0, -m (1, 1) / d);
      res.set (1, 1, -1.0 / d);
      break;
    case 'S': // H to S
      res = htos (m);
      break;
    }
    break;
  case 'G':
    switch (out) {
    case 'Y': // G to Y
      d = m (1, 1);
      res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
      res.set (0, 1, m (0, 1) / d);
      res.set (1, 0, -m (1, 0) / d);
      res.set (1, 1, 1.0 / d);
      break;
    case 'Z': // G to Z
      d = m (0, 0);
      res.set (0, 0, 1.0 / d);
      res.set (0, 1, -m (0, 1) / d);
      res.set (1, 0, m (1, 0) / d);
      res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d);
      break;
    case 'H': // G to H
      d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
      res.set (0, 0, m (1, 1) / d);
      res.set (0, 1, -m (0, 1) / d);
      res.set (1, 0, -m (1, 0) / d);
      res.set (1, 1, m (0, 0) / d);
      break;
    case 'G': // G to G
      res = m;
      break;
    case 'A': // G to A
      d = m (1, 0);
      res.set (0, 0, 1.0 / d);
      res.set (0, 1, m (1, 1) / d);
      res.set (1, 0, m (0, 0) / d);
      res.set (1, 1, m (0, 0) * m (1, 1) / d - m (0, 1));
      break;
    case 'S': // G to S
      res = gtos (m);
      break;
    }
    break;
  case 'A':
    switch (out) {
    case 'Y': // A to Y
      d = m (0, 1);
      res.set (0, 0, m (1, 1) / d);
      res.set (0, 1, m (1, 0) - m (0, 0) * m (1, 1) / d);
      res.set (1, 0, -1.0 / d);
      res.set (1, 1, m (0, 0) / d);
      break;
    case 'Z': // A to Z
      d = m (1, 0);
      res.set (0, 0, m (0, 0) / d);
      res.set (0, 1, m (0, 0) * m (1, 1) / d - m (0, 1));
      res.set (1, 0, 1.0 / d);
      res.set (1, 1, m (1, 1) / d);
      break;
    case 'H': // A to H
      d = m (1, 1);
      res.set (0, 0, m (0, 1) / d);
      res.set (0, 1, m (0, 0) - m (0, 1) * m (1, 0) / d);
      res.set (1, 0, -1.0 / d);
      res.set (1, 1, m (1, 0) / d);
      break;
    case 'G': // A to G
      d = m (0, 0);
      res.set (0, 0, m (1, 0) / d);
      res.set (0, 1, m (1, 0) * m (0, 1) / d - m (1, 1));
      res.set (1, 0, 1.0 / d);
      res.set (1, 1, m (0, 1) / d);
      break;
    case 'A': // A to A
      res = m;
      break;
    case 'S': // A to S
      res = atos (m);
      break;
    }
    break;
  case 'S':
    switch (out) {
    case 'S': // S to S
      res = m;
      break;
    case 'T': // S to T
      d = m (1, 0);
      res.set (0, 0, m (0, 1) - m (0, 0) * m (1, 1) / d);
      res.set (0, 1, m (0, 0) / d);
      res.set (1, 0, -m (1, 1) / d);
      res.set (0, 1, 1.0 / d);
      break;
    case 'A': // S to A
      res = stoa (m);
      break;
    case 'H': // S to H
      res = stoh (m);
      break;
    case 'G': // S to G
      res = stog (m);
      break;
    case 'Y': // S to Y
      res = stoy (m);
      break;
    case 'Z': // S to Z
      res = stoz (m);
      break;
    }
    break;
  case 'T':
    switch (out) {
    case 'S': // T to S
      d = m (1, 1);
      res.set (0, 0, m (0, 1) / d);
      res.set (0, 1, m (0, 0) - m (0, 1) * m (1, 0) / d);
      res.set (1, 0, 1.0 / d);
      res.set (0, 1, -m (1, 0) / d);
      break;
    case 'T': // T to T
      res = m;
      break;
    }
    break;
  }
  return res;
}

nr_double_t rollet (matrix m) {
  assert (m.getRows () >= 2 && m.getCols () >= 2);
  nr_double_t res;
  res = (1 - norm (m (0, 0)) - norm (m (1, 1)) + norm (det (m))) /
    2 / abs (m (0, 1) * m (1, 0));
  return res;
}

/* Computes stability measure B1 of the given S-parameter matrix. */
nr_double_t b1 (matrix m) {
  assert (m.getRows () >= 2 && m.getCols () >= 2);
  nr_double_t res;
  res = 1 + norm (m (0, 0)) - norm (m (1, 1)) - norm (det (m));
  return res;
}

} // namespace qucs