---

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

---