matrix.cc

模糊聚类分析的源程序！

      for (i = na; i >= 0; i--) {
	w = Lese_i_j (i, i) - p;
	r = Lese_i_j (i, en);
	for (j = m; j <= na; j++)
	  r += Lese_i_j (i, j) * Lese_i_j (j, en);
	if (wi[i] < 0.0) {
	  z = w;
	  s = r;
	} else {
	  m = i;
	  if (wi[i] == 0.0) {
	    if (w != 0.0)
	      Setze_i_j (i, en, -r / w);
	    else
	      Setze_i_j (i, en, -r / (MACH_EPS * norm));
	  } else {
	    x = Lese_i_j (i, i + 1);
	    y = Lese_i_j (i + 1, i);
	    q = wr[i] - p;
	    q = q * q + wi[i] * wi[i];
	    t = (x * s - z * r) / q;
	    Setze_i_j (i, en, t);
	    if (fabs (x) > fabs (z))
	      Setze_i_j (i + 1, en, (-r - w * t) / x);
	    else
	      Setze_i_j (i + 1, en, (-s - y * t) / z);
	  }
	}
      }
    } else {
      /* komplexe Loesung */
      Fehlermeldung ("DMatrix::hqrvec(int low, int high, DVektor& wr, DVektor& wi, DMatrix& eivec))", KOMPLEXFEHLER);
    }
  }

  for (i = 0; i < n; i++) {
    if ((i < low) || (i > high))
      for (k = i + 1; k < n; k++)
	eivec.Setze_i_j (i, k, Lese_i_j (i, k));
  }
  for (j = n - 1; j >= low; j--) {
    if (j < high)
      m = j;
    else
      m = high;
    if (wi[j] < 0.0) {
      l = j - 1;
      for (i = low; i <= high; i++) {
	y = 0.0;
	z = 0.0;
	for (k = low; k <= m; k++) {
	  y += eivec.Lese_i_j (i, k) * Lese_i_j (k, l);
	  z += eivec.Lese_i_j (i, k) * Lese_i_j (k, j);
	}
	eivec.Setze_i_j (i, l, y);
	eivec.Setze_i_j (i, j, z);
      }
    } else {
      if (wi[j] == 0.0) {
	for (i = low; i <= high; i++) {
	  z = 0.0;
	  for (k = low; k <= m; k++)
	    z += eivec[i][k] * Lese_i_j (k, j);
	  eivec.Setze_i_j (i, j, z);
	}
      }
    }
  }
  return 0;
}

char 
DMatrix::Hqr2 (int low, int high, DVektor & wr,
	       DVektor & wi, DMatrix & eivec, IVektor & cnt)
{
  int na, en, iter, i, j, k, l, m;
  int enditer, stop;
  double p = 0, q = 0, r = 0, s, t, w, x, y, z;

  for (i = 0; i < n; i++) {
    if ((i < low) || (i > high)) {
      wr[i] = Lese_i_j (i, i);
      wi[i] = 0.0;
      cnt[i] = 0;
    }
  }
  en = high;
  t = 0.0;

  while (en >= low) {
    iter = 0;
    na = en - 1;
    enditer = FALSE;
    do {
      l = en;
      while (l > low) {
	if (fabs (Lese_i_j (l, l - 1)) <=
	MACH_EPS * (fabs (Lese_i_j (l - 1, l - 1)) + fabs (Lese_i_j (l, l))))
	  goto cont;
	l--;
      }

    cont:x = Lese_i_j (en, en);
      if (l == en) {
	Setze_i_j (en, en, x + t);
	wr[en] = Lese_i_j (en, en);
	wi[en] = 0.0;
	cnt[en] = iter;
	en--;
	enditer = TRUE;
      } else {
	y = Lese_i_j (na, na);
	w = Lese_i_j (en, na) * Lese_i_j (na, en);
	if (l == na) {
	  p = (y - x) * 0.5;
	  q = p * p + w;
	  z = sqrt (fabs (q));
	  x = x + t;
	  Setze_i_j (en, en, x);
	  Setze_i_j (na, na, y + t);
	  cnt[en] = -iter;
	  cnt[na] = iter;
	  if ((q < 0) && (q + 1E-5 > 0))
	    q = 0;
	  if (q >= 0.0) {
	    if (p < 0.0)
	      z = p - z;
	    else
	      z = p + z;
	    wr[na] = x + z;
	    s = x - w / z;
	    wr[en] = s;
	    wi[na] = 0.0;
	    wi[en] = 0.0;
	    x = Lese_i_j (en, na);
	    r = sqrt (x * x + z * z);
	    p = x / r;
	    q = z / r;
	    for (j = na; j < n; j++) {
	      z = Lese_i_j (na, j);
	      Setze_i_j (na, j, q * z + p * Lese_i_j (en, j));
	      Setze_i_j (en, j, q * Lese_i_j (en, j) - p * z);
	    }
	    for (i = 0; i <= en; i++) {
	      z = Lese_i_j (i, na);
	      Setze_i_j (i, na, q * z + p * Lese_i_j (i, en));
	      Setze_i_j (i, en, q * Lese_i_j (i, en) - p * z);
	    }
	    for (i = low; i <= high; i++) {
	      z = eivec[i][na];
	      eivec[i][na] = q * z + p * eivec[i][en];
	      eivec[i][en] = q * eivec[i][en] - p * z;
	    }
	  }
	  /* q>=0.0 */ 
	  else {
	    Fehlermeldung ("DMatrix::Hqr2(int low, int high,DVektor& wr, DVektor& wi, DMatrix& eivec, IVektor& cnt)", KOMPLEXFEHLER);
	    wr[na] = x + p;
	    wr[en] = x + p;
	    wi[na] = z;
	    wi[en] = -z;
	  }
	  en = en - 2;
	  enditer = TRUE;
	} else {		/* l==na */
	  enditer = FALSE;
	  if (iter >= MAXIT) {
	    cnt[en] = 1 + MAXIT;
	    return en;
	  }
	  if ((iter != 0) && ((iter % KORITER) == 0)) {
	    t += x;
	    for (i = low; i <= en; i++)
	      Setze_i_j (i, i, Lese_i_j (i, i) - x);
	    s = fabs (Lese_i_j (en, na)) + fabs (Lese_i_j (na, en - 2));
	    y = 0.75 * s;
	    x = y;
	    w = -0.4375 * s * s;
	  }
	  iter++;
	  m = en - 1;
	  stop = (en - 2 < l);
	  while ((m >= l) && (!stop)) {
	    m--;
	    z = Lese_i_j (m, m);
	    r = x - z;
	    s = y - z;
	    p = (r * s - w) / Lese_i_j (m + 1, m) + Lese_i_j (m, m + 1);
	    q = Lese_i_j (m + 1, m + 1) - z - r - s;
	    r = Lese_i_j (m + 2, m + 1);
	    s = fabs (p) + fabs (q) + fabs (r);
	    p = p / s;
	    q = q / s;
	    r = r / s;
	    stop = (m <= l);
	    if (m > l)
	      stop = (fabs (Lese_i_j (m, m - 1)) * (fabs (q) + fabs (r)) <=
		      MACH_EPS * fabs (p) * (fabs (Lese_i_j (m - 1, m - 1)) + fabs (z) + fabs (Lese_i_j (m + 1, m + 1))));
	  }
	  for (i = m + 2; i <= en; i++)
	    Setze_i_j (i, i - 2, 0.0);
	  for (i = m + 3; i <= en; i++)
	    Setze_i_j (i, i - 3, 0.0);
	  for (k = m; k <= na; k++) {
	    if (k != m) {
	      p = Lese_i_j (k, k - 1);
	      q = Lese_i_j (k + 1, k - 1);
	      if (k != na)
		r = Lese_i_j (k + 2, k - 1);
	      else
		r = 0.0;
	      x = fabs (p) + fabs (q) + fabs (r);
	      if (x == 0.0)
		goto weiter;
	      p = p / x;
	      q = q / x;
	      r = r / x;
	    }
	    s = sqrt (p * p + q * q + r * r);
	    if (p < 0.0)
	      s = -s;
	    if (k != m)
	      Setze_i_j (k, k - 1, -s * x);
	    else if (l != m)
	      Setze_i_j (k, k - 1, -Lese_i_j (k, k - 1));
	    p = p + s;
	    x = p / s;
	    y = q / s;
	    z = r / s;
	    q = q / p;
	    r = r / p;
	    for (j = k; j < n; j++) {
	      p = Lese_i_j (k, j) + q * Lese_i_j (k + 1, j);
	      if (k != na) {
		p = p + r * Lese_i_j (k + 2, j);
		Setze_i_j (k + 2, j, Lese_i_j (k + 2, j) - p * z);
	      }
	      Setze_i_j (k + 1, j, Lese_i_j (k + 1, j) - p * y);
	      Setze_i_j (k, j, Lese_i_j (k, j) - p * x);
	    }
	    if (k + 3 < en)
	      j = k + 3;
	    else
	      j = en;
	    for (i = 0; i <= j; i++) {
	      p = x * Lese_i_j (i, k) + y * Lese_i_j (i, k + 1);
	      if (k != na) {
		p += z * Lese_i_j (i, k + 2);
		Setze_i_j (i, k + 2, Lese_i_j (i, k + 2) - p * r);
	      }
	      Setze_i_j (i, k + 1, Lese_i_j (i, k + 1) - p * q);
	      Setze_i_j (i, k, Lese_i_j (i, k) - p);
	    }
	    for (i = low; i <= high; i++) {
	      p = x * eivec.Lese_i_j (i, k) + y * eivec.Lese_i_j (i, k + 1);
	      if (k != na) {
		p = p + z * eivec.Lese_i_j (i, k + 2);
		eivec.Setze_i_j (i, k + 2, eivec.Lese_i_j (i, k + 2) - p * r);
	      }
	      eivec.Setze_i_j (i, k + 1, eivec.Lese_i_j (i, k + 1) - p * q);
	      eivec.Setze_i_j (i, k, eivec.Lese_i_j (i, k) - p);
	    }
	  weiter:;
	  }			/* k */
	}			/* l!=na */
      }				/* l!=en */
    } while (!enditer);
  }
  if (Hqrvec (low, high, wr, wi, eivec) != 0)
    return 99;
  return 0;
};

void 
DMatrix::Elmtrans (int low, int high, IVektor & perm, DMatrix & Result)
{
  int i, j, k;

  Result.Einheitsmatrix ();
  for (i = high - 1; i >= low + 1; i--) {
    j = perm[i];
    for (k = i + 1; k <= high; k++)
      Result.Setze_i_j (k, i, Lese_i_j (k, i - 1));
    if (i != j) {
      for (k = i; k <= high; k++) {
	Result.Setze_i_j (i, k, Result.Lese_i_j (j, k));
	Result.Setze_i_j (j, k, 0.0);
      }
      Result.Setze_i_j (j, i, 1.0);
    }
  }
}

char 
DMatrix::Elmhes (int low, int high, IVektor & perm)
{
  int i, j, m;
  double x, y;

  for (m = low + 1; m < high; m++) {
    i = m;
    x = 0.0;
    for (j = m; j <= high; j++) {
      if (fabs (Lese_i_j (j, m - 1)) > fabs (x)) {
	x = Lese_i_j (j, m - 1);
	i = j;
      }
    }

    perm[m] = i;
    if (i != m) {
      for (j = m - 1; j < n; j++)
	swap (Lese_i_j (i, j), Lese_i_j (m, j));
      for (j = 0; j <= high; j++)
	swap (Lese_i_j (j, i), Lese_i_j (j, m));
    }
    if (x != 0.0) {
      for (i = m + 1; i <= high; i++) {
	y = Lese_i_j (i, m - 1);
	if (y != 0.0) {
	  Setze_i_j (i, m - 1, y / x);
	  y = Lese_i_j (i, m - 1);
	  for (j = m; j < n; j++)
	    Lese_i_j (i, j) -= y * Lese_i_j (m, j);
	  for (j = 0; j <= high; j++)
	    Lese_i_j (j, m) += y * Lese_i_j (j, i);
	}
      }
    }
  }
};

char 
DMatrix::Balance (DVektor & Skal, int &low, int &high)
{
  int i, j, k, m;
  char contin;
  double b2, r, c, f, g, s;

  if (this->n != this->m) {
    Fehlermeldung ("DMatrix::Balance(DVektor& Skal, int& low, int& high)", DIMFEHLER);
    return (1);
  }
  Skal.Setze_Dim (this->n);
  b2 = BASIS * BASIS;
  m = 0;
  k = n - 1;
  do {
    contin = FALSE;
    for (j = k; j >= 0; j--) {
      r = 0.0;
      for (i = 0; i <= k; i++) {
	if (i != j)
	  r += fabs (Lese_i_j (j, i));
      }
      if (r == 0.0) {
	Skal[k] = j;
	if (j != k) {
	  for (i = 0; i <= k; i++)
	    swap (Lese_i_j (i, j), Lese_i_j (i, k));
	  for (i = m; i < n; i++)
	    swap (Lese_i_j (j, i), Lese_i_j (k, i));
	}
	k--;
	contin = TRUE;
      }
    }
  } while (contin);

  do {
    contin = FALSE;
    for (j = m; j <= k; j++) {
      c = 0.0;
      for (i = m; i <= k; i++) {
	if (i != j)
	  c += fabs (Lese_i_j (i, j));
      }
      if (c == 0.0) {
	Skal[m] = j;
	if (j != m) {
	  for (i = 0; i <= k; i++)
	    swap (Lese_i_j (i, j), Lese_i_j (i, m));
	  for (i = m; i < n; i++)
	    swap (Lese_i_j (j, i), Lese_i_j (m, i));
	}
	m++;
	contin = TRUE;
      }
    }
  } while (contin);

  low = m;
  high = k;
  for (i = m; i <= k; i++)
    Skal[i] = 1.0;
  do {
    contin = FALSE;
    for (i = m; i <= k; i++) {
      c = 0.0;
      r = 0.0;
      for (j = m; j <= k; j++)
	if (j != i) {
	  c += fabs (Lese_i_j (j, i));
	  r += fabs (Lese_i_j (i, j));
	}
      g = r / BASIS;
      f = 1.0;
      s = c + r;
      while (c < g) {
	f *= BASIS;
	c *= b2;
      }
      g = r * BASIS;
      while (c >= g) {
	f /= BASIS;
	c /= b2;
      }
      if (((c + r) / f) < (0.95 * s)) {
	g = 1.0 / f;
	Skal[i] = Skal[i] * f;
	contin = TRUE;
	for (j = m; j < n; j++)
	  Lese_i_j (i, j) *= g;
	for (j = 0; j <= k; j++)
	  Lese_i_j (j, i) *= f;
      }
    }
  } while (contin);

  return (0);
};

void 
DMatrix::Speichern (FILE * File)
{
  int i, j, k = 0;

  fprintf (File, "%d %d\n", m, n);
  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++)
      fprintf (File, "%lf ", Daten[k++]);
    fprintf (File, "\n");
  }
};

int 
DMatrix::Laden (FILE * File)
{
  int i, j, k = 0;

  if (!(fscanf (File, "%d %d", &i, &j))) {
#if 0
    this. ~ DMatrix ();
#endif
    return (0);
  }
  Setze_Dim (i, j);
  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++)
      if (!(fscanf (File, "%lf", &Daten[k++]))) {
#if 0
	this. ~ DMatrix ();
#endif
	return (0);
      }
  }
  return (1);
};

void 
DMatrix::Quicksort (int a, int b, int Index)
{
  int r, i, j;
  double x, DBuffer;

  if (b > a) {
    i = a;
    j = b;
    /* Waehle r einmal in der Mitte */
    r = (a + b) / 2;
    x = Lese_i_j (r, Index);
    do {
      while (Lese_i_j (i, Index) < x)
	i++;
      while (Lese_i_j (j, Index) > x)
	j--;
      if (i <= j) {
	/* tauschen */