qavrgmatrix.cpp

Go to the documentation of this file.

#include "qavrgmatrix.h"

#include <math.h>
#include <stdio.h>

int QavrgMatrix::m_Errors = 0;

QavrgMatrix::QavrgMatrix(int nr, int nc, QObject *parent)
  : QObject(parent),
  m_Matrix(nr*nc),
  m_NRows(nr),
  m_NCols(nc)
{
}

void QavrgMatrix::resize(int nr, int nc)
{
  m_NRows = nr;
  m_NCols = nc;

  m_Matrix.resize(nr*nc);
}

double& QavrgMatrix::operator()(int r, int c)
{
  if (r >= 0 && r < m_NRows && c >= 0 && c < m_NCols) {
    return m_Matrix[r*m_NCols+c];
  } else {
    printf("Array access (%d,%d) outside bounds (%d,%d)\n", r,c, m_NRows, m_NCols);
    return m_Dummy;
  }
}

int QavrgMatrix::size(int d)
{
  if (d==1) return m_NRows;
  if (d==2) return m_NCols;
  return 0;
}

QavrgMatrix& QavrgMatrix::operator=(const QavrgMatrix& m)
                                   {
  m_Matrix = m.m_Matrix;

  return *this;
}

void QavrgMatrix::fill(double val)
{
  m_Matrix.fill(val);
}

void QavrgMatrix::dump()
{
  printf("QavrgMatrix = (%d,%d) {\n", m_NRows, m_NCols);

  for (int i=0; i<m_NRows; i++) {
    if (i) printf(",\n");
    printf(" { ");
    for (int j=0; j<m_NCols; j++) {
      if (j) printf(",\t");
      printf("%g",operator()(i,j));
    }
    printf(" }");
  }

  printf("\n}\n");
}

void QavrgMatrix::test()
{
  QavrgMatrix m(5,5);

  for (int i=0; i<5; i++) {
    for (int j=0; j<5; j++) {
      m(i,j) = 10*i+j;
    }
  }

  m.dump();
}

#define SWAP(a,b) { double temp=(a); (a)=(b); (b)=temp; }

void QavrgMatrix::gaussj(QavrgMatrix &a, QVector<double> &b)
{
  int n1= a.size(1);
  int n2= a.size(2);
  int n = b.size();

  if (n == n1 && n == n2) {
    int irow, icol;
    QVector<int> indxc(n), indxr(n), ipiv(n);
    ipiv.fill(-1);

    for (int i=0; i<n; i++) {
      double big=0;

      for (int j=0; j<n; j++) {
        if (ipiv[j] != 0) {
          for (int k=0; k<n; k++) {
            if (ipiv[k] == -1) {
              if (fabs(a(j,k)) >= big) {
                big = fabs(a(j,k));
                irow=j;
                icol=k;
              }
            } else if (ipiv[k] > 0) {
              if (m_Errors < 10) {
                m_Errors++;
                printf("gaussj: Singular Matrix-1\n");
              }

              return;
            }
          }
        }
      }

      ++(ipiv[icol]);
      
      if (irow != icol) {
        for (int l=0; l<n; l++) {
          SWAP(a(irow,l),a(icol,l));
        }
        SWAP(b[irow],b[icol]);
      }

      indxr[i]=irow;
      indxc[i]=icol;

      if (a(icol,icol) == 0.0) {
        if (m_Errors < 10) {
          m_Errors++;
          printf("gaussj: Singular Matrix-2\n");
        }

        return;
      }

      double pivinv = 1.0/a(icol,icol);
      a(icol,icol) = 1.0;

      for (int l=0; l<n; l++) {
        a(icol,l) *= pivinv;
      }
      b[icol] *= pivinv;

      for (int ll=0; ll<n; ll++) {
        if (ll != icol) {
          double dum = a(ll,icol);
          a(ll,icol) = 0.0;

          for (int l=0; l<n; l++) {
            a(ll,l) -= a(icol,l)*dum;
          }

          b[ll] -= b[icol]*dum;
        }
      }
    }

    for (int l=n-1; l>=0; l--) {
      if (indxr[l] != indxc[l]) {
        for (int k=0; k<n; k++) {
          SWAP(a(k,indxr[l]), a(k,indxc[l]));
        }
      }
    }
  }
}