---

matrix_invert.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 inverts a (non-singular) matrix.                 */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
MATRIX_INVERT (double *mat,     /*  I   Pointer to the input matrix[size][size].              */
               short size,      /*  I   Size of the matrix.                                   */
               double *inv      /*  O   Pointer to the inverted matrix[size][size].           */
/*----------------------------------------------------------------------------*/
  )
{
  register short i, j, k;
  double temp;

  if (TIMES)
    TIMING (T_MATRIX_INVERT);

  /* copy imput matrix into output matrix */
  for (i = 0; i < size; i++)
    {
      for (j = 0; j < size; j++)
        {
          inv[i * size + j] = mat[i * size + j];
        }
    }

  /* perform inversion on inv, in place */
  for (i = 0; i < size; i++)
    {
      temp = inv[i * size + i];
      inv[i * size + i] = 1.0;
      for (j = 0; j < size; j++)
        {
          if (temp == 0.0)
            {
              MESSAGE ('E', " Can't invert singular matrix.");
              goto the_end;
            }
          else
            {
              inv[i * size + j] /= temp;
            }
        }
      for (j = 0; j < size; j++)
        {
          if (j != i)
            {
              temp = inv[j * size + i];
              inv[j * size + i] = 0.0;
              for (k = 0; k < size; k++)
                {
                  inv[j * size + k] -= temp * inv[i * size + k];
                }
            }
        }
    }

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

---