---

eigen.c

Go to the documentation of this file.

#include        "sadie.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This low-level function computes the eigenvalues and eigenvectors of a   */
/*   symmetric matrix.                                                        */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Background Information-----------------------------------------------------*/
/*                                                                            */
/*   Davis, J.C.:                                                             */
/*   "Statistics and Data Analysis in Geology."                               */
/*   J. Wiley & Sons, New York, 1986                                          */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
EIGEN (double *mat,             /*  I   Pointer to the symmetric input matrix[size][size].    */
       short size,              /*  I   Size of the matrix.                                   */
       double *eig              /*  O   Pointer to the eigenvector matrix[size][size].        */
/*----------------------------------------------------------------------------*/
  )
{
  register short i, j, k, found;
  char msg[SLEN];
  double norm, thr, al, am, ao, sine, sinsq, cosine, cossq, at, bt, xt, temp;

  if (TIMES)
    TIMING (T_EIGEN);

  /* set eig to identity matrix */
  for (norm = 0.0, i = 0; i < size; i++)
    {
      for (j = 0; j < size; j++)
        {
          if (i == j)
            {
              eig[i * size + j] = 1.0;
            }
          else
            {
              eig[i * size + j] = 0.0;
              norm += mat[i * size + j] * mat[i * size + j];
            }
        }
    }

  /* iterate */
  thr = sqrt (norm);
  norm = 1.0E-09 * sqrt (norm) / (double) size;
  do
    {
      thr /= (double) size;
      do
        {

          /* scan down columns for off-diagonal elements greater than or equal to threshold */
          for (found = FALSE, i = 1; i < size; i++)
            {
              for (j = 0; j < i; j++)
                {
                  if (thr <= fabs (mat[j * size + i]))
                    {

                      /* compute sines and cosines */
                      found = TRUE;
                      al = -mat[j * size + i];
                      am = (mat[j * size + j] - mat[i * size + i]) / 2.0;
                      ao =
                        am <
                        0.0 ? -al / sqrt (al * al +
                                          am * am) : al / sqrt (al * al +
                                                                am * am);
                      sine = ao / sqrt (2.0 * (1.0 + sqrt (1.0 - ao * ao)));
                      sinsq = sine * sine;
                      cosine = sqrt (1.0 - sinsq);
                      cossq = cosine * cosine;

                      /* rotate columns i and j */
                      for (k = 0; k < size; k++)
                        {
                          if (k != i && k != j)
                            {
                              at = mat[k * size + j];
                              mat[k * size + j] =
                                at * cosine - mat[k * size + i] * sine;
                              mat[k * size + i] =
                                at * sine + mat[k * size + i] * cosine;
                            }
                          bt = eig[k * size + j];
                          eig[k * size + j] =
                            bt * cosine - eig[k * size + i] * sine;
                          eig[k * size + i] =
                            bt * sine + eig[k * size + i] * cosine;
                        }
                      at = mat[j * size + j];
                      bt = mat[i * size + i];
                      xt = 2.0 * mat[j * size + i] * sine * cosine;
                      mat[j * size + j] = at * cossq + bt * sinsq - xt;
                      mat[i * size + i] = at * sinsq + bt * cossq + xt;
                      mat[i * size + j] = mat[j * size + i] =
                        (at - bt) * sine * cosine + mat[j * size +
                                                        i] * (cossq - sinsq);
                      for (k = 0; k < size; k++)
                        {
                          mat[j * size + k] = mat[k * size + j];
                          mat[i * size + k] = mat[k * size + i];
                        }
                    }
                }
            }
        }
      while (found);
    }
  while (thr > norm);

  /* sort eigenvalues and eigenvectors */
  for (i = 1; i < size; i++)
    {
      for (j = i; j > 0; j--)
        {
          if (mat[(j - 1) * size + j - 1] < mat[j * size + j])
            {
              temp = mat[(j - 1) * size + j - 1];
              mat[(j - 1) * size + j - 1] = mat[j * size + j];
              mat[j * size + j] = temp;
              for (k = 0; k < size; k++)
                {
                  temp = eig[k * size + j - 1];
                  eig[k * size + j - 1] = eig[k * size + j];
                  eig[k * size + j] = temp;
                }
            }
        }
    }

  /* output eigenvalues and eigenvectors */
  MESSAGE ('I', "");
  MESSAGE ('I', " Eigenvalues        Eigenvectors");
  for (i = 0; i < size; i++)
    {
      sprintf (msg, " %d   %11.4e     ", i + 1, mat[i * size + i]);
      for (j = 0; j < size; j += 1)
        {
          sprintf (msg + strlen (msg), "%12.4e", eig[j * size + i]);
        }
      MESSAGE ('I', msg);
    }

the_end:
  if (TIMES)
    TIMING (T_EXIT);
}

---