matrix.cc

模糊聚类分析的源程序！

/********************************************************/
/*   filename: matrix.cc                                */
/*                                                      */
/********************************************************/
/* programmed by: Oliver Wagner                         */
/* last change: 26-04-95                                */
/********************************************************/
#include "matrix.h"

/********************************************************/
/*              Alles zur Klasse DMatrix                */
/********************************************************/
DMatrix::DMatrix (int mGroesse, int nGroesse, double *Inhalt)
{
  int i;

  n = nGroesse;
  m = mGroesse;
  if ((n != 0) && (m != 0)) {
    if ((Daten = (double *) malloc (sizeof (double) * m * n)) == NULL)
        Fehlermeldung ("DMatrix", SPEICHERFEHLER);

    if (Inhalt != NULL)
      for (i = 0; i < n * m; i++)
	Daten[i] = Inhalt[i];
    else
      for (i = 0; i < n * m; i++)
	Daten[i] = 0;
  } else
    Daten = NULL;
};

DMatrix::DMatrix (const DMatrix & Daraus)
{
  int i;

  n = Daraus.n;
  m = Daraus.m;
  if ((n != 0) && (m != 0)) {
    if ((Daten = (double *) malloc (sizeof (double) * m * n)) == NULL)
        Fehlermeldung ("DMatrix::DMatrix(const DMatrix& Daraus)", SPEICHERFEHLER);

    for (i = 0; i < n * m; i++)
      Daten[i] = Daraus.Daten[i];
  } else
    Daten = NULL;
};

DMatrix::~DMatrix ()
{
  if ((n != 0) && (m != 0)) {
    n = m = 0;
    free (Daten);
    Daten = NULL;
  }
};

void 
DMatrix::print ()
{
  int i, j;

  printf ("\nDMatrix :");
  for (i = 0; i < m; i++) {
    if (i != 0)
      printf ("         ");	/* Versatz */
    for (j = 0; j < n; j++)
      printf ("%lf ", Daten[i * n + j]);
    printf ("\n");
  }
};

void 
DMatrix::Setze_Dim (int Neues_m, int Neues_n)
{

  if ((m != Neues_m) || (n != Neues_n)) {
    if ((m != 0) && (n != 0))
      free (Daten);
    m = Neues_m;
    n = Neues_n;
    if ((Daten = (double *) malloc (sizeof (double) * m * n)) == NULL)
        Fehlermeldung ("DMatrix", SPEICHERFEHLER);
  }
};

void 
DMatrix::Einheitsmatrix (int Groesse)
{
  int i;

  Setze_Dim (Groesse, Groesse);

  for (i = 0; i < Groesse * Groesse; i++)
    Daten[i] = 0;
  for (i = 0; i < Groesse * Groesse; i += Groesse + 1)
    Daten[i] = 1;
};

void 
DMatrix::Einheitsmatrix ()
{
  int i;

  for (i = 0; i < m * n; i++)
    Daten[i] = 0;
  for (i = 0; i < m * n; i += n + 1)
    Daten[i] = 1;
};

void 
DMatrix::NullMatrix ()
{
  int i;

  for (i = 0; i < n * m; i++)
    Daten[i] = 0;
};

DMatrix 
DMatrix::Projiziere (BVektor & Diese)
return Result (Diese.Anzahl_gesetzt (), Diese.Anzahl_gesetzt (), NULL);
{
  int i, j, a, b;

  a = 0;
  for (i = 0; i < m; i++) {
    b = 0;
    if (Diese[i]) {
      for (j = 0; j < n; j++) {
	if (Diese[j])
	  Result[a][b++] = (*this)[i][j];
      }
      a++;
    }
  }
}

DMatrix 
DMatrix::Projiziere (IVektor & Diese)
return Result (Diese.Lese_Dim (), Diese.Lese_Dim (), NULL);
{
  int i, j;

  for (i = 0; i < Result.m; i++) {
    for (j = 0; j < Result.n; j++) {
      Result[i][j] = (*this)[Diese[i]][Diese[j]];
    }
  }
}

DMatrix 
DMatrix::Deprojiziere (BVektor & Diese)
return Result (Diese.Lese_Dim (), Diese.Lese_Dim (), NULL);
{
  int i, j, a, b;

  a = 0;
  for (i = 0; i < Result.m; i++) {
    if (Diese[i]) {
      b = 0;
      for (j = 0; j < Result.n; j++) {
	if (Diese[j])
	  Result[i][j] = (*this)[a][b++];
      }
      a++;
    }
  }
}

DMatrix 
DMatrix::Gauss_Inverse ()
return Result (m, m, NULL);
{
  DMatrix *Buffer;
  int Spalte, Zeile;
  double Faktor;
  #define Inverse_eps 1e-12
  #define Inverse2_eps 1e-50


  Result.Einheitsmatrix (m);
  Buffer = new DMatrix (*this);

  for (Spalte = 0; Spalte < m; Spalte++) {
    if ((Buffer->Lese_i_j (Spalte, Spalte) <= Inverse_eps) &&
	(Buffer->Lese_i_j (Spalte, Spalte) >= -Inverse_eps)) {	/* moeliche Zeile zum Tauschen suchen */
      for (Zeile = Spalte + 1; Zeile < m; Zeile++) {
	if ((Buffer->Lese_i_j (Zeile, Spalte) >= Inverse_eps) ||
	    (Buffer->Lese_i_j (Zeile, Spalte) <= -Inverse_eps)) {
	  Buffer->Tausche_Zeilen (Spalte, Zeile);
	  Result.Tausche_Zeilen (Spalte, Zeile);
	  break;
	}
      }

      /* nullwert in letzter zeile/spalte */
      if (Zeile == m && 
	  (Buffer->Lese_i_j (Spalte, Spalte) <= Inverse2_eps) &&
	  (Buffer->Lese_i_j (Spalte, Spalte) >= -Inverse2_eps)) {
	printf("%g\n",Buffer->Lese_i_j (Spalte, Spalte));
	Fehlermeldung ("DMatrix::Gauss_Inverse()", INVERTIERFEHLER);
	return (Result);
      } else if (Zeile==m) {
	printf("%g\n",Buffer->Lese_i_j (Spalte, Spalte));
	printf ("Warning: matrix inversion might be unreliable\n");
      }
    }
    for (Zeile = 0; Zeile < m; Zeile++) {
      if ((Zeile != Spalte) &&
	  ((Buffer->Lese_i_j (Zeile, Spalte) >= Inverse_eps) ||
	   (Buffer->Lese_i_j (Zeile, Spalte) <= -Inverse_eps))) {
	/* Aufraeumen */
	Faktor = -Buffer->Lese_i_j (Zeile, Spalte) /
	  Buffer->Lese_i_j (Spalte, Spalte);

	Buffer->Zeile_plus_Zeile (Zeile, Spalte, Faktor);
	Result.Zeile_plus_Zeile (Zeile, Spalte, Faktor);
      }
    }
  }
  /* Diagonale auf 1 bringen */
  for (Spalte = 0; Spalte < m; Spalte++)
    Result.Zeile_Multiplizieren (Spalte, 1 / Buffer->Lese_i_j (Spalte, Spalte));

  delete Buffer;
};

double 
DMatrix::Determinante ()
{
  double Result = 1;
  DMatrix *Buffer;
  int Spalte, Zeile;
  double Faktor;

  Buffer = new DMatrix (*this);

  for (Spalte = 0; Spalte < m; Spalte++) {
    if (Buffer->Lese_i_j (Spalte, Spalte) == 0) {
      Result = 0;
      break;
    } else {
      for (Zeile = Spalte + 1; Zeile < m; Zeile++) {
	/* Aufraeumen */
	Faktor = -Buffer->Lese_i_j (Zeile, Spalte) /
	  Buffer->Lese_i_j (Spalte, Spalte);
	Buffer->Zeile_plus_Zeile (Zeile, Spalte, Faktor);
      }
    }
    /* Diagonale aufmultiplizieren */
    for (Spalte = 0; Spalte < m; Spalte++)
      Result *= Buffer->Lese_i_j (Spalte, Spalte);
  }

  delete Buffer;
  return (Result);
};

double 
DMatrix::Zeilensumme (int i)
{
  int x, bis;
  double Result = 0;

  bis = x + n;
  for (x = i * n; x < (i + 1) * n; x++)
    Result += Daten[x];
  return (Result);
};

void 
DMatrix::Tausche_Zeilen (int x, int y)
{
  double Buffer;
  int i;

  for (i = 0; i < n; i++) {
    Buffer = Lese_i_j (x, i);
    Setze_i_j (x, i, Lese_i_j (y, i));
    Setze_i_j (y, i, Buffer);
  }
};

void 
DMatrix::Tausche_Spalten (int x, int y)
{
  double Buffer;
  int i;

  for (i = 0; i < m; i++) {
    Buffer = Lese_i_j (i, x);
    Setze_i_j (i, x, Lese_i_j (i, y));
    Setze_i_j (i, y, Buffer);
  }
};

void 
DMatrix::Zeile_Multiplizieren (int Zeile, double Damit)
{
  int i;

  for (i = Zeile * n; i < (Zeile + 1) * n; i++)
    Daten[i] *= Damit;
};

void 
DMatrix::Spalte_Multiplizieren (int Spalte, double Damit)
{
  int i;

  for (i = Spalte; i < n * m; i += n)
    Daten[i] *= Damit;
};

void 
DMatrix::Zeile_plus_Zeile (int Diese, int plusDiese, double Faktor)
{
  int i;

  for (i = 0; i < n; i++)
    Setze_i_j (Diese, i, Lese_i_j (Diese, i) + Lese_i_j (plusDiese, i) * Faktor);
};

DVektor 
DMatrix::Spalte_aus_Matrix (int i)
  return Result (m, i, NULL);
{
  int j;

  for (j = 0; j < m; j++)
    Result[j] = Lese_i_j (j, i);
};

DVektor 
DMatrix::Zeile_aus_Matrix (int i)
  return Result (n, i, NULL);
{
  int j;
  for (j = 0; j < n; j++)
    Result[j] = Lese_i_j (i, j);
};

void 
DMatrix::Setze_Spalte (int Die, DVektor & Damit)
{
  int i;

  for (i = 0; i < m; i++)
    Daten[Die + i * n] = Damit[i];
};

DVektorArray 
DMatrix::Eigenvektoren ()
{
  IVektor Iterationen (m, NULL);

  Iterationen = ITERATIONEN;
  return (Eigenvektoren (Iterationen));
}

DVektorArray 
DMatrix::Eigenvektoren (IVektor & cnt)
return Result (n, n, NULL);
{
  DVektor Skal, valreal (n, 0, NULL), valim (n, 0, NULL);
  DMatrix Uebergabe_Matrix (*this), Eivec (n, n, NULL);
  int low, high, i;

  if (n < 2)
    Fehlermeldung ("DMatrix::Eigenvektoren(IVektor& cnt)", DIMFEHLER);

  Uebergabe_Matrix.Balance (Skal, low, high);
  Uebergabe_Matrix.Elmhes (low, high, cnt);
  Uebergabe_Matrix.Elmtrans (low, high, cnt, Eivec);
  Uebergabe_Matrix.Hqr2 (low, high, valreal, valim,
			 Eivec, cnt);
  Eivec.Balback (low, high, Skal);
  Eivec.Norm_1 (valim);

  for (i = 0; i < n; i++)
    Result[i] = Eivec.Spalte_aus_Matrix (i);
};

DMatrix 
DMatrix::Eigenmatrix ()
{
  IVektor Iterationen (m, NULL);

  Iterationen = ITERATIONEN;
  return (Eigenmatrix (Iterationen));
}

DMatrix 
DMatrix::Eigenmatrix (IVektor & cnt)
return Result (n, n, NULL);
{
  DVektor Skal, valreal (n, 0, NULL), valim (n, 0, NULL);
  DMatrix Uebergabe_Matrix (*this);
  int low, high;

  if (n < 2)
    Fehlermeldung ("DMatrix::Eigenvektoren(IVektor& cnt)", DIMFEHLER);

  Uebergabe_Matrix.Balance (Skal, low, high);
  Uebergabe_Matrix.Elmhes (low, high, cnt);
  Uebergabe_Matrix.Elmtrans (low, high, cnt, Result);
  Uebergabe_Matrix.Hqr2 (low, high, valreal, valim,
			 Result, cnt);
  Result.Balback (low, high, Skal);
  Result.Norm_1 (valim);
};

void 
DMatrix::Norm_1 (DVektor & wi)
{
  int i, j;
  double maxi;

  j = 0;
  while (j < n) {
    if (wi[j] == 0.0) {
      maxi = Lese_i_j (0, j);
      for (i = 1; i < n; i++) {
	if (fabs (Lese_i_j (i, j)) > fabs (maxi))
	  maxi = Lese_i_j (i, j);
      }
      if (maxi != 0.0) {
	maxi = 1.0 / maxi;
	for (i = 0; i < n; i++)
	  Lese_i_j (i, j) *= maxi;
      }
      j++;
    } else {
      Fehlermeldung ("DMatrix::Norm_1(DMatrix wi)", KOMPLEXFEHLER);
    }
  }
}

void 
DMatrix::Balback (int low, int high, DVektor & skal)
{
  int i, j, k;
  double s;

  for (i = low; i <= high; i++) {
    s = skal[i];
    for (j = 0; j < n; j++)
      Lese_i_j (i, j) *= s;
  }

  for (i = low - 1; i >= 0; i--) {
    k = (int) floor (skal[i] + 0.4999);
    if (k != i)
      Tausche_Zeilen (i, k);
  }

  for (i = high + 1; i < n; i++) {
    k = (int) floor (skal[i] + 0.4999);
    if (k != i)
      Tausche_Zeilen (i, k);
  }
}

char 
DMatrix::Hqrvec (int low, int high, DVektor & wr, DVektor & wi,
		 DMatrix & eivec)
{
  int na, en, i, j, k, l, m;
  double norm, p, q, r, s = 0, t, w, x, y, z = 1;

  norm = 0.0;
  k = 0;
  for (i = 0; i < n; i++) {
    for (j = k; j < n; j++)
      norm += fabs (Lese_i_j (i, j));
    k = i;
  }
  if (norm == 0.0)
    return 1;			/* Nullmatrix */

  for (en = n - 1; en >= 0; en--) {
    p = wr[en];
    q = wi[en];
    na = en - 1;
    if (q == 0.0) {
      m = en;
      Setze_i_j (en, en, 1.0);