home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATCHOL.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  3.8 KB  |  129 lines
  1. /**************************************************/
  2. /*     matchol.c source for Cholesky Functions    */
  3. /**************************************************/
  4.  
  5.  
  6. /**************************************************/
  7. /*            MatClass Source File                */
  8. /*       Copyright of C. R. Birchenhall           */
  9. /*       University of Manchester, UK.            */
  10. /*   MatClass is freeware. This file should be    */
  11. /* made freely available to users of any software */
  12. /* whose creation is wholly or partly dependent   */
  13. /*                on this file.                   */
  14. /**************************************************/
  15.  
  16. #include "matchol.hpp"
  17.  
  18. /***************************************************/
  19. /*                   cholesky                      */
  20. /***************************************************/
  21.  
  22. #include <math.h>
  23. // for sqrt
  24.  
  25. INDEX cholesky( matrix& x, matError& error  )
  26. /********************************************
  27.   Cholesky decomposition. Lower triang of x
  28.   overwritten by decomposition. Nonsquare
  29.   and non-definite matrices trapped as
  30.   errors.
  31. ********************************************/
  32. {
  33.    static char *mName = "cholesky" ;
  34.    matFunc func( mName ) ; x.debugInfo( func ) ;
  35.    INDEX i, j, k, n = x.nRows() ;
  36.    DOUBLE g, h ;
  37.    error = NOERROR ;
  38.    if ( n != x.nCols() ) {
  39.       error = NTSQR ;
  40.       return 0 ;
  41.    } // if
  42.    for ( j = 1 ; j <= n; j++ ) {
  43.       g = x(j,j) ;
  44.       for ( k = 1 ; k < j ; k++ )
  45.         g -= x(j,k) * x(j,k) ;
  46.       if ( g <= 0.0 ) {
  47.          error = NOTPD ;
  48.          return 0 ;
  49.       } // if
  50.       x(j,j) = g = sqrt( g ) ;
  51.       for ( i = j + 1 ; i <= n ; i++ ) {
  52.          h = x(i,j) ;
  53.          for ( k = 1 ; k < j ; k++ )
  54.             h -= x(i,k) * x(j,k) ;
  55.          x(i,j) = h / g ;
  56.       } // for i
  57.    } // for j
  58.    return OK ;
  59. } // cholesky
  60.  
  61. INDEX cholSolve( const matrix& A, matrix& b, REAL tol,
  62.                  matError& error )
  63. /****************************************************************
  64.   Solve equations Ax=b assuming A has been factored by
  65.   cholesky. The lower triang factor G is assumed to be
  66.   in lower triang of A. The real tol is used to determine
  67.   if diagonal element is zero. The solution is returned in b.
  68. *****************************************************************/
  69. {
  70.    static char *mName = "cholSolve" ;
  71.    matFunc func( mName ) ; A.debugInfo( func ) ;
  72.    INDEX i, j, n = A.nRows() ;
  73.    DOUBLE r ;
  74.    error = NOERROR ;
  75.    if ( n != A.nCols() ) {
  76.       error = NTSQR ;
  77.       return 0 ;
  78.    } // if
  79.    if ( n != b.nRows() ) {
  80.       error = NEDIM ;
  81.       return 0 ;
  82.    } // if
  83.    // Use forward substitution on G
  84.    for ( i = 1 ; i <= n ; i++ ) {
  85.       if ( A(i,i) < tol ) {
  86.          error = NOTPD ;
  87.          return 0 ;
  88.       } // if
  89.       r = b(i) ;
  90.       for ( j = 1 ; j < i ; j++ )
  91.          r -= A(i,j) * b(j) ;
  92.       b(i) = r / A(i,i) ;
  93.    } // for i
  94.    // Use backward substitution on G transpose
  95.    for ( i = n ; i ; i-- ) {
  96.       // no need to check diagonals again
  97.       r = b(i) ;
  98.       for ( j = n ; j > i ; j-- )
  99.          r -= A(j,i) * b(j) ;
  100.       b(i) = r / A(i,i) ;
  101.    } // for i
  102.    return OK ;
  103. } // cholSolve
  104.  
  105. matrix normalEqn( const matrix& x, const matrix& b )
  106. {
  107.    static char *mName = "normalEqn" ;
  108.    matFunc func( mName ) ; x.debugInfo( func ) ;
  109.    if ( b.nRows() != x.nRows() )
  110.       b.errorExit( mName, NEDIM ) ;
  111.    matError error ;
  112.    REAL tol = matrixTol() ;
  113.    INDEX nr = x.nCols(), nc = b.nCols() ;
  114.    matrix z( nr, nc ) ;
  115.    z.TMultOf( x, b ) ;
  116.    matrix xTx( nr, nr ) ; ;
  117.    xTx.TMultOf( x, x ) ;
  118.    cholesky( xTx, error ) ;
  119.    refMatrix zc( z ) ;
  120.    for ( INDEX i = 1; !error && i <= nc ; i++ ) {
  121.       zc.refCol(i) ;
  122.       cholSolve( xTx, zc, tol, error ) ;
  123.    } //for
  124.    if ( error )
  125.       xTx.errorExit( mName, error ) ;
  126.    return z ;
  127. } // normalEqn
  128.  
  129.