---

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);
}

---