Power-Programmierung/CD2.mdf/c/library/dos/diverses/leda/src/basic/_matrix.c

C/C++ Source or Header | 1991-11-15 | 7.1 KB | 360 lines

/*******************************************************************************
+
+  LEDA  2.1.1                                                 11-15-1991
+
+
+  _matrix.c
+
+
+  Copyright (c) 1991  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 6600 Saarbruecken, FRG
+  All rights reserved.
+
*******************************************************************************/



/*  Vektoren und Matrizen */

#include <LEDA/matrix.h>

#include <math.h>

#define ELEM(M,i,j)  M.ptr->v[i].ptr->v[j]

#define EPSILON 1e-12

//------------------------------------------------------------------------------
//  matrix member functions
//------------------------------------------------------------------------------


matrix_rep::matrix_rep(int d1, int d2)
{
 if (d1<0 || d2<0)
  error_handler(1,"matrix: negative dimension.");

  dim1=d1;
  dim2=d2;

  if (dim1 > 0)
  { v = (vector*) new char*[dim1];
    while (d1--) v[d1].ptr = new vector_rep(d2);
   }
  else v = nil;
}

matrix_rep::matrix_rep(const matrix_rep* p)
{
  int d1 = (p) ? p->dim1 : 0;
  int d2 = (p) ? p->dim2 : 0;

  dim1=d1;
  dim2=d2;

  if (dim1 > 0)
  { v = (vector*) new char*[dim1];
    while (d1--) v[d1].ptr = new vector_rep(p->v[d1].ptr);
   }
  else v = nil;
}

matrix_rep::matrix_rep(int d1, int d2, double** p)
{
  dim1=d1;
  dim2=d2;
  v = (vector*) new char*[dim1];
  while (d1--)
  { v[d1].ptr = new vector_rep(dim2);
    for(int i=0;i<d2;i++) v[d1].ptr->v[i] = p[d1][i];
   }

 }

matrix_rep::~matrix_rep()
{ if (v)
  { while (dim1--) delete v[dim1].ptr;
    delete (char*)v;
   }
 }

void matrix::check_dimensions(const matrix& mat) const
{ if (ptr->dim1 != mat.ptr->dim1 || ptr->dim2 != mat.ptr->dim2)
   error_handler(1,"incompatible matrix types.");
 }

matrix::matrix(const vector& vec)
{ int d1 = vec.ptr->dim;
  ptr = new matrix_rep(d1,1);
  while (d1--) ptr->v[d1].ptr->v[0] = vec.ptr->v[d1];
}

matrix& matrix::operator=(const matrix& mat)
{ register int d1 = mat.ptr->dim1;
  register int d2 = mat.ptr->dim2;
  if (ptr->dim1 != d1 || ptr->dim2 != d2)
  { clear();
    ptr = new matrix_rep(d1,d2);
   }
  while (d1--) ptr->v[d1] = mat.ptr->v[d1];
  return *this;
}

int matrix::operator==(const matrix& x) const
{ int d1 = x.ptr->dim1;
  int d2 = x.ptr->dim2;

  if (ptr->dim1 != d1 || ptr->dim2 != d2) return false;

  for(int i=0;i<d1;i++)
    for(int j=0;j<d2;j++)
      if (elem(i,j) != x.elem(i,j)) return false;

  return true;
}


vector& matrix::row(int i) const
{ if ( i<0 || i>=ptr->dim1 )  error_handler(1,"matrix: row index out of range");
   return ptr->v[i];
}


double& matrix::operator()(int i, int j)
{ if ( i<0 || i>=ptr->dim1 )  error_handler(1,"matrix: row index out of range");
  if ( j<0 || j>=ptr->dim2 )  error_handler(1,"matrix: col index out of range");
  return ptr->v[i].ptr->v[j];
}

double  matrix::operator()(int i, int j) const
{ if ( i<0 || i>=ptr->dim1 )  error_handler(1,"matrix: row index out of range");
  if ( j<0 || j>=ptr->dim2 )  error_handler(1,"matrix: col index out of range");
  return ptr->v[i].ptr->v[j];
}

vector matrix::col(int i)  const
{ if ( i<0 || i>=ptr->dim2 )  error_handler(1,"matrix: col index out of range");
  vector result(ptr->dim1);
  int j = ptr->dim1;
  while (j--) result.ptr->v[j] = elem(j,i);
  return result;
}

matrix::operator vector() const
{ if (ptr->dim2!=1)
   error_handler(1,"error: cannot make vector from matrix\n");
  return col(0);
}

matrix matrix::operator+(const matrix& mat)
{ check_dimensions(mat);

  matrix result(ptr->dim1,ptr->dim2);
  register int n = ptr->dim1;
  while (n--) result.ptr->v[n] = ptr->v[n]+mat.ptr->v[n];
  return result;
}

matrix matrix::operator-(const matrix& mat)
{ check_dimensions(mat);

  matrix result(ptr->dim1,ptr->dim2);
  register int n = ptr->dim1;
  while (n--) result.ptr->v[n] = ptr->v[n]-mat.ptr->v[n];
  return result;
}

matrix matrix::operator*(double f)
{ matrix result(ptr->dim1,ptr->dim2);
  register int n = ptr->dim1;
  while (n--) result.ptr->v[n] = ptr->v[n]*f;
  return result;
}

matrix matrix::operator*(const matrix& mat)
{ if (ptr->dim2!=mat.ptr->dim1)
     error_handler(1,"matrix multiplication: incompatible matrix types\n");

  matrix result(ptr->dim1, mat.ptr->dim2);
  register int i;
  register int j;

  for (i=0;i<mat.ptr->dim2;i++)
  for (j=0;j<ptr->dim1;j++) ELEM(result,j,i) = ptr->v[j] * mat.col(i);

 return result;

}

double matrix::det() const
{
 if (ptr->dim1!=ptr->dim2)
   error_handler(1,"matrix::det: matrix not quadratic.\n");

 int n = ptr->dim1;

 matrix M(n,1);

 int flips;

 double** A = triang(M,flips);

 double Det = 1;

 for(int i=0;i<n;i++) Det *= A[i][i];

 for(i=0;i++;i<n) delete A[i];
 delete A;

 return (flips % 2) ? -Det : Det;

}


double** matrix::triang(const matrix& M, int& flips)  const
{
 register double **p, **q;
 register double *l, *r, *s;

 register double pivot_el,tmp;

 register int i,j, col, row;

 int n = ptr->dim1;
 int d = M.ptr->dim2;
 int m = n+d;

 double** A = new double*[n];

 p = A;

 for(i=0;i<n;i++)
 { *p = new double[m];
   l = *p++;
   for(j=0;j<n;j++) *l++ = elem(i,j);
   for(j=0;j<d;j++) *l++ = M.elem(i,j);
  }

 flips = 0;

 for (col=0, row=0; row<n; row++, col++)
 {
   // search for row j with maximal entry in current col
   j = row;
   for (i=row+1; i<n; i++)
     if (A[j][col] < A[i][col]) j = i;

   if ( n > j && j > row)
   { SWAP(A[j],A[row]);
     flips++;
    }

   tmp = A[row][col];
   q  = &A[row];

   if (fabs(tmp) < EPSILON) error_handler(1,"matrix has not full rank.");

   for (p = &A[n-1]; p != q; p--)
   {
     l = *p+col;
     s = *p+m;
     r = *q+col;

     if (*l != 0.0)
     { pivot_el = *l/tmp;
        while(l < s) *l++ -= *r++ * pivot_el;
      }

    }

  }

 return A;
}

matrix matrix::inv() const
{
 if (ptr->dim1!=ptr->dim2)
     error_handler(1,"matrix::inv: matrix not quadratic.\n");
 int n = ptr->dim1;
 matrix I(n,n);
 for(int i=0; i<n; i++) I(i,i) = 1;
 return solve(I);
}



matrix matrix::solve(const matrix& M) const
{

if (ptr->dim1 != ptr->dim2 || ptr->dim1 != M.ptr->dim1)
     error_handler(1, "Solve: wrong dimensions\n");

 register double **p, ** q;
 register double *l, *r, *s;

 int      n = ptr->dim1;
 int      d = M.ptr->dim2;
 int      m = n+d;
 int      row, col,i;


 double** A = triang(M,i);

 for (col = n-1, p = &A[n-1]; col>=0; p--, col--)
 {
   s = *p+m;

   double tmp = (*p)[col];

   for(l=*p+n; l < s; l++) *l /=tmp;

   for(q = A; q != p; q++ )
   { tmp = (*q)[col];
     l = *q+n;
     r = *p+n;
     while(r < s)  *l++ -= *r++ * tmp;
    }

  }

  matrix result(n,d);

  for(row=0; row<n; row++)
  { l = A[row]+n;
    for(col=0; col<d; col++) ELEM(result,row,col) = *l++;
    delete A[row];
   }

  delete A;

  return result;
}




matrix matrix::trans() const
{ matrix result(ptr->dim2,ptr->dim1);
  for(int i = 0; i < ptr->dim2; i++)
    for(int j = 0; j < ptr->dim1; j++)
      ELEM(result,i,j) = elem(j,i);
  return result;
}


ostream& operator<<(ostream& s, const matrix& M)
{ int i;
  s <<"\n";
  for (i=0;i<M.dim1();i++) s << M[i] << "\n";
  return s;
}

istream& operator>>(istream& s, matrix& M)
{ int i=0;
  while (i<M.dim1() && s >> M[i++]);
  return s;
}