00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include <stdio.h>
00008 #include <math.h>
00009 #include <stdlib.h>
00010 #include "eispack.h"
00011 
00012 /*--------------------------------------------------------------------------
00013    Compute the eigenvalue/vector decomposition of a symmetric matrix,
00014    stored in double precision.
00015      n = order of matrix
00016      a = on input: matrix(i,j) is in a[i+n*j] for i=0..n-1 , j=0..n-1
00017            output: a[i+n*j] has the i'th component of the j'th eigenvector
00018      e = on input: not used (but the calling program must
00019                              allocate the space for e[0..n-1])
00020            output: e[j] has the j'th eigenvalue, ordered so that
00021            e[0] <= e[1] <= ... <= e[n-1]
00022    Uses the f2c translation of EISPACK.
00023 ----------------------------------------------------------------------------*/
00024 
00025 void symeig_double( int n , double *a , double *e )
00026 {
00027    integer nm , matz , ierr ;
00028    double *fv1 , *fv2 ;
00029 
00030    if( a == NULL || e == NULL || n < 1 ) return ;
00031 
00032    if( n == 1 ){
00033      e[0] = a[0] ; a[0] = 1.0 ; return ;  /* degenerate case */
00034    }
00035 
00036    fv1 = (double *) malloc(sizeof(double)*n) ;  /* workspaces */
00037    fv2 = (double *) malloc(sizeof(double)*n) ;
00038 
00039    nm = n ; matz = 1 ; ierr = 0 ;
00040 
00041    rs_( &nm , &nm , a , e , &matz , a , fv1 , fv2 , &ierr ) ;
00042 
00043    free((void *)fv1) ; free((void *)fv2) ;
00044    return ;
00045 }
00046 
00047 /*------------------ just compute the eigenvalues -------------------*/
00048 
00049 void symeigval_double( int n , double *a , double *e )
00050 {
00051    integer nm , matz , ierr ;
00052    double *fv1 , *fv2 ;
00053 
00054    if( a == NULL || e == NULL || n < 1 ) return ;
00055 
00056    if( n == 1 ){ e[0] = a[0] ; return ; } /* degenerate case */
00057 
00058    fv1 = (double *) malloc(sizeof(double)*n) ;  /* workspaces */
00059    fv2 = (double *) malloc(sizeof(double)*n) ;
00060 
00061    nm = n ; matz = 0 ; ierr = 0 ;
00062 
00063    rs_( &nm , &nm , a , e , &matz , a , fv1 , fv2 , &ierr ) ;
00064 
00065    free((void *)fv1) ; free((void *)fv2) ;
00066    return ;
00067 }
00068 
00069 /*--------------------------------------------------------------------*/
00070 /*! Compute SVD of double precision matrix:                      T
00071                                             [a] = [u] diag[s] [v]
00072     - m = # of rows in a
00073     - n = # of columns in a
00074     - a = pointer to input matrix; a[i+j*m] has the (i,j) element
00075     - s = pointer to output singular values; length = n
00076     - u = pointer to output matrix, if desired; length = m*n
00077     - v = pointer to output matrix, if desired; length = n*n
00078 
00079 ----------------------------------------------------------------------*/
00080 
00081 void svd_double( int m , int n , double *a , double *s , double *u , double *v )
00082 {
00083    integer mm,nn , lda,ldu,ldv , ierr ;
00084    doublereal *aa, *ww , *uu , *vv , *rv1 ;
00085    logical    matu , matv ;
00086 
00087    if( a == NULL || s == NULL || m < 1 || n < 1 ) return ;
00088 
00089    mm  = m ;
00090    nn  = n ;
00091    aa  = a ;
00092    lda = m ;
00093    ww  = s ;
00094 
00095    if( u == NULL ){
00096      matu = (logical) 0 ;
00097      uu   = (doublereal *)malloc(sizeof(double)*m*n) ;
00098    } else {
00099      matu = (logical) 1 ;
00100      uu = u ;
00101    }
00102    ldu = m ;
00103 
00104    if( v == NULL ){
00105      matv = (logical) 0 ;
00106      vv   = NULL ;
00107    } else {
00108      matv = (logical) 1 ;
00109      vv = v ;
00110    }
00111    ldv = n ;
00112 
00113    rv1 = (double *) malloc(sizeof(double)*n) ;  /* workspace */
00114 
00115    (void) svd_( &mm , &nn , &lda , aa , ww ,
00116                 &matu , &ldu , uu , &matv , &ldv , vv , &ierr , rv1 ) ;
00117 
00118    free((void *)rv1) ;
00119 
00120    if( u == NULL ) free((void *)uu) ;
00121    return ;
00122 }
00123 
00124 #if 0
00125 /*--------------------------------------------------------------------*/
00126 /*! Compute SVD of single precision matrix:                      T
00127                                             [a] = [u] diag[s] [v]
00128     - m = # of rows in a
00129     - n = # of columns in a
00130     - a = pointer to input matrix; a[i+j*m] has the (i,j) element
00131     - s = pointer to output singular values; length = n
00132     - u = pointer to output matrix, if desired; length = m*n
00133     - v = pointer to output matrix, if desired; length = n*n
00134 
00135 ----------------------------------------------------------------------*/
00136 
00137 void svd_float( int m , int n , float *a , float *s , float *u , float *v )
00138 {
00139    integer mm,nn , lda,ldu,ldv , ierr ;
00140    doublereal *aa, *ww , *uu , *vv , *rv1 ;
00141    logical    matu , matv ;
00142 
00143    int ii ;
00144 
00145    if( a == NULL || s == NULL || m < 1 || n < 1 ) return ;
00146 
00147    mm  = m ;
00148    nn  = n ;
00149    aa  = (doublereal *)malloc(sizeof(double)*m*n) ;
00150    lda = m ;
00151    ww  = (doublereal *)malloc(sizeof(double)*n) ;
00152 
00153    for( ii=0 ; ii < m*n ; ii++ ) aa[ii] = (doublereal)a[ii] ;
00154 
00155    matu = (logical) 1 ;
00156    uu   = (doublereal *)malloc(sizeof(double)*m*n) ;
00157    ldu  = m ;
00158 
00159    matv = (logical) 1 ;
00160    vv   = (doublereal *)malloc(sizeof(double)*n*n) ;
00161    ldv  = n ;
00162 
00163    rv1 = (double *) malloc(sizeof(double)*n) ;  /* workspace */
00164 
00165    (void) svd_( &mm , &nn , &lda , aa , ww ,
00166                 &matu , &ldu , uu , &matv , &ldv , vv , &ierr , rv1 ) ;
00167 
00168    free((void *)rv1) ; free((void *)aa) ;
00169 
00170    /* copy results out */
00171 
00172    for( ii=0 ; ii < n ; ii++ ) s[ii] = (float)ww[ii] ;
00173 
00174    if( u != NULL )
00175      for( ii=0 ; ii < m*n ; ii++ ) u[ii] = (float)uu[ii] ;
00176 
00177    if( v != NULL )
00178      for( ii=0 ; ii < n*n ; ii++ ) v[ii] = (float)vv[ii] ;
00179 
00180    free((void *)uu) ; free((void *)vv) ;
00181    return ;
00182 }
00183 #endif