 /*
  * Khoros: $Id: matrix.c,v 1.2 1992/03/20 23:42:51 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: matrix.c,v 1.2 1992/03/20 23:42:51 dkhoros Exp $";
#endif

 /*
  * $Log: matrix.c,v $
 * Revision 1.2  1992/03/20  23:42:51  dkhoros
 * VirtualPatch5
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * 
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *----------------------------------------------------------------------
 */

/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   >>>>                                                       <<<<
   >>>>	    file name: matrix.c                               <<<<
   >>>>                                                       <<<<
   >>>>   description: Contains routines that are used for    <<<<
   >>>>                simple matrix and vector operations.   <<<<
   >>>>                                                       <<<<
   >>>>                This is part of vmath.a                <<<<
   >>>>                                                       <<<<
   >>>>      routines: 					      <<<<
   >>>>                                                       <<<<
   >>>> modifications: See module headers                     <<<<
   >>>>                                                       <<<<
   >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */


#include "unmcopyright.h"	 /* Copyright 1991 by UNM */
#include "vinclude.h"	

/***********************************************************************
*
*  Routine Name: INNER
*
*          Date: 4-Jan-1991
*
*       Purpose: Compute inner product of two vectors
*
*         Input: x - (double *) to first vector 
*                y - (double *) to second vector
*                n - number of components in the vectors
*
*        Output: (double)x'y
*
*    Written By:  Scott Wilson
*
* Modifications: 
*
***********************************************************************/
double
inner(x,y,n)
double *x,*y;
int n;
  {
    int i;
    double a;

    a = 0;
    for (i=0; i<n; i++) a += x[i]*y[i];
    return(a);
  }

/***********************************************************************
*
*  Routine Name: MATVEC
*
*          Date: 4-Jan-1991
*
*       Purpose: Compute matrix-vector product
*
*         Input: a - (double *) to matrix 
*                x - (double *) to input vector
*                y - (double *) to output vector
*                rows - number of rows in matrix a
*                cols - number of cols in matrix a
*
*        Output: (double)y=Ax
*
*    Written By: Scott Wilson
*
* Modifications: 
*
***********************************************************************/
matvec(y,a,x,rows,cols)
double *a,*x,*y;
int rows,cols;
  {
    int i;
    double inner();

    for (i=0; i<rows; i++) y[i] = inner(a+i*cols,x,cols);
  }

/***********************************************************************
*
*  Routine Name: GET_PRINC_AXIS
*
*          Date: 4-Jan-1991
*
*       Purpose: Obtain the principal axis of a covariance matrix
*
*         Input: a - (double *) to matrix 
* 		 n - size of matrix (nXn)
*                y - (double *) to output vector (axis)
*
*        Output: (double)y eigenvector associated with the dominant
*		 eigenvalue of a.
*
*    Written By:  Scott Wilson
*
* Modifications: 
*
***********************************************************************/
/*
    get_princ_axis() - Obtain the principal axis (the eigenvector
    associated with the largest eigenvalue) for the covariance matrix
    "a" of size nXn. Put the principal eigenvector in the place pointed
    to by "y". NOTE: a and y should be of type DOUBLE! Also note that
    since a is a covariance matrix it is symmetric, real, and very
    likely to be positive definate. It is positive semidefinate for sure,
    and may also be diagonal dominant.

    How this thing works:
      Use the power iteration method to obtain the dominant eigenvalue
      and its associated eigenvector as described in Gollub and VanLoan,
      MATRIX COMPUTATIONS, pp 209.

    Possible difficulties: The major hitch in this technique is that the
    convergence is proportional to lambda(2)/lambda(1) where lambda(1)
    is the dominant eigenvector. If these eigenvalues are closely spaced
    then we won't get a decent eigenvector (it will have an incorrect
    direction).

    Fortunately, when using the principal axis to split a cluster and there
    are two very strongly dominant axes  with the same ellipticity, then
    we can split on any combination of those axis and reduce the cluster
    variances greatly.

    The 10 iterations used in the code have been found to be satisfactory
    for all of the data I have encountered unless I fix up a really nasty
    special case.
*/
get_princ_axis(a,n,y)
double *a,*y;
int n;
  {
    double *b,*c,inner(),lambda;
    int i,j;
    char *malloc();

    c = (double *)malloc(n*sizeof(double));
    b = (double *)malloc(n*sizeof(double));
    
    for (i=0; i<n; i++) c[i] = 1.0;
    for (j=0; j<10; j++)
      {
        matvec(b,a,c,n,n);
        lambda = sqrt(inner(b,b,n)); 
        for (i=0; i<n; i++) c[i] = b[i]/lambda;
      }
    for (i=0; i<n; i++) y[i] = c[i];
  }

