home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATLU.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  10.4 KB  |  349 lines
  1. /**************************************************/
  2. /*         matlu.c source for LU 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 "matlu.hpp"
  17. #include <math.h>
  18.  
  19. /****************************************************
  20.             crout.cxx
  21.  
  22.  Crout's Method of LU decomposition of square matrix,
  23.  with implicit partial pivoting.
  24.  A is overwritten : U is explicit in the upper triag
  25.  and L is in multiplier form in the subdiagionals
  26.  i.e. subdiag A.mat(i,j) is the multiplier used to
  27.  eliminate the (i,j) term. Row permutations are
  28.  mapped out in index. sign is used for calculating
  29.  the determinant, see determinant function below.
  30.  
  31.  Adapted from Numerical Recipes in C by Press et al.
  32. *****************************************************/
  33.  
  34. INDEX crout( matrix& A, indexArray& index, REAL& sign, REAL tol,
  35.              matError& error  )
  36. {
  37.    static char *mName = "crout" ;
  38.    matFunc func( mName ) ; A.debugInfo( func ) ;
  39.    INDEX n = A.nRows() ;
  40.    error = NOERROR ;
  41.    if ( n != A.nCols() ) {
  42.       error = NTSQR ;
  43.       return 0 ;
  44.    } // if
  45.    if ( n != index.length() )
  46.       index.reset( n ) ;
  47.  
  48.    INDEX    i, imax, j, k, zeros = 0 ;
  49.    DOUBLE   sum ;
  50.    REAL     big, temp ;
  51.    matrix   scale(n) ;
  52.  
  53.    sign = 1.0 ;
  54.  
  55.    /*****************************
  56.     Find implicit scaling factors
  57.    ******************************/
  58.    for ( i = 1; i <= n ; i++ ) {
  59.       big = 0.0 ;
  60.       for ( j = 1; j <= n; j++ ) {
  61.          temp = (REAL) fabs( (double) A.mat(i,j) ) ;
  62.          if ( temp > big )
  63.             big = temp ;
  64.       } // for
  65.       scale(i) = ( big == 0.0 ? 0.0 : 1.0 / big ) ;
  66.    } // for i
  67.  
  68.    for ( j = 1; j <= n; j++ ) {
  69.       /*************************************
  70.        Run down jth column from top to diag,
  71.        to form the elements of U.
  72.       **************************************/
  73.       for ( i = 1; i < j; i++ ) {
  74.          sum = (DOUBLE) A(i,j) ;
  75.          for ( k = 1; k < i; k++ )
  76.             sum -= A(i,k) * A(k,j) ;
  77.          A(i,j) = (REAL) sum ;
  78.       } // for i
  79.       /******************************************
  80.        Run down jth subdiag to form the residuals
  81.        after the elimination of the first j-1
  82.        subdiags.  These residuals divided by the
  83.        appropriate diagonal term will become the
  84.        multipliers in the elimination of the jth.
  85.        subdiag. Find index of largest scaled term
  86.        in imax.
  87.       *******************************************/
  88.       big = 0.0 ;
  89.       for ( i = j; i <= n; i++ ) {
  90.          sum = (DOUBLE) A(i,j) ;
  91.          for ( k = 1; k < j; k++ )
  92.             sum -= A(i,k) * A(k,j) ;
  93.          A(i,j) = (REAL) sum ;
  94.          temp = scale(i) * ( (REAL) fabs( sum ) ) ;
  95.          if ( temp >= big ) {
  96.             big = temp ;
  97.             imax = i ;
  98.          } // if
  99.       }  // for i
  100.       /*****************************
  101.        Permute current row with imax
  102.       ******************************/
  103.       if ( j != imax ) {
  104.          for ( k = 1; k <= n; k++ ) {
  105.             temp = A( imax, k ) ;
  106.             A( imax, k ) = A(j,k) ;
  107.             A(j,k) = temp ;
  108.          } // for k
  109.          sign = - sign ;
  110.          scale( imax ) = scale(j) ;
  111.       } // if
  112.       index(j) = imax;
  113.       /***************************************
  114.        If diag term is not zero divide subdiag
  115.        to form multipliers.
  116.       ****************************************/
  117.       if ( fabs( A(j,j) ) < tol )
  118.          zeros++ ;
  119.       else if ( j != n ) {
  120.          temp = 1.0 / A(j,j) ;
  121.          for ( i = j+1; i <= n; i++ )
  122.             A(i,j) *= temp;
  123.       } // if
  124.    } // for j
  125.    if ( zeros )
  126.       error = SINGM ;
  127.    return zeros ;
  128. } // crout()
  129.  
  130.  
  131. /***************************************************
  132.                     luSolve.cxx
  133.  
  134.  Solve Ax=b assuming A is in the LU form, but
  135.  assume b has *not* been transformed.  Solution
  136.  returned in b.  A is unchanged.
  137.  
  138.  Adapted from Numerical Recipes in C by Press et al.
  139. ****************************************************/
  140.  
  141. INDEX luSolve( const matrix& A, const indexArray& index,
  142.            matrix& b, REAL tol, matError& error )
  143. {
  144.    static char *mName = "luSolve" ;
  145.    matFunc func( mName ) ; A.debugInfo( func ) ;
  146.    INDEX  n = A.nRows() ;
  147.    error = NOERROR ;
  148.    if ( b.nCols() != 1 ) {
  149.       error = NOTVC ;
  150.       return 0 ;
  151.    } // if
  152.    if ( n != b.nRows() ) {
  153.       error = NEDIM ;
  154.       return 0 ;
  155.    } // if
  156.  
  157.    INDEX  i, nonzero=0, iperm, j ;
  158.    DOUBLE sum, zero = 0.0 ;
  159.  
  160.    /*************************
  161.     Check for zero diagonals
  162.    **************************/
  163.    for ( i = 1; i <= n ; i++ )
  164.       if ( fabs( A(i,i) ) < tol )
  165.          return i ;
  166.  
  167.    /**************************************
  168.     Transform b allowing for leading zeros
  169.    ***************************************/
  170.  
  171.    for ( i = 1; i <= n; i++ ) {
  172.       iperm = index(i) ;
  173.       sum = (DOUBLE) b( iperm ) ;
  174.       b( iperm ) = b(i) ;
  175.       if ( nonzero )
  176.          for ( j = nonzero; j <= i-1; j++ )
  177.             sum -= A(i,j) * b(j) ;
  178.       else if ( sum != zero )
  179.          nonzero = i ;
  180.       b(i) = (REAL) sum ;
  181.    } // for i
  182.  
  183.    /****************
  184.     backsubstitution
  185.    *****************/
  186.  
  187.    for ( i = n; i >= 1; i-- ) {
  188.       sum = (DOUBLE) b(i) ;
  189.       for ( j = i+1; j <= n; j++ )
  190.          sum -= A(i,j) * b(j) ;
  191.       b(i) = (REAL) ( sum / A(i,i) ) ;
  192.    } // for i
  193.    return 0 ;
  194. } // luSolve
  195.  
  196.  
  197.  
  198. /***************************************************
  199.                   luInverse.cxx
  200.  
  201.  Returns inverse of A in invA assuming A in LU form.
  202. ****************************************************/
  203.  
  204.  
  205. INDEX luInverse( const matrix& A, const indexArray& index,
  206.          matrix& invA, REAL tol, matError& error )
  207. {
  208.    static char *mName = "luInverse" ;
  209.    matFunc func( mName ) ; A.debugInfo( func ) ;
  210.    INDEX  nc = A.nCols(), nr = A.nRows() ;
  211.    if ( ( nc != invA.nCols() ) || ( nr != invA.nRows() ) )
  212.       invA.reset( nr, nc ) ;
  213.    INDEX j, exit = 0 ;
  214.    invA = 0.0 ;
  215.    invA.setDiag( 1.0 ) ;
  216.    refMatrix col( invA ) ;
  217.    error = NOERROR ;
  218.    for ( j = 1 ; ( exit == 0 ) && !error && ( j <= nc ) ; j++ ) {
  219.       col.refCol(j) ;
  220.       exit = luSolve( A, index, col, tol, error ) ;
  221.    } // for
  222.    return exit ;
  223. } // luInverse
  224.  
  225.  
  226. /***************************************************
  227.                     lutSolve.cxx
  228.  
  229.  Solve A.trans()x=b assuming A is in the LU form,
  230.  but assume b has *not* been transformed.
  231.  Solution returned in b.  A is unchanged.
  232.  
  233.  Adapted from luSolve above and based on LINPACK
  234.  USER's GUIDE p1.14.
  235. ****************************************************/
  236.  
  237. INDEX lutSolve( const matrix& A, const indexArray& index,
  238.             matrix& b, REAL tol, matError& error )
  239. {
  240.    static char *mName = "lutSolve" ;
  241.    matFunc func( mName ) ; A.debugInfo( func ) ;
  242.    INDEX  n = A.nCols() ;
  243.    error = NOERROR ;
  244.    if ( b.nCols() != 1 ) {
  245.       error = NOTVC ;
  246.       return 0 ;
  247.    } // if
  248.    if ( n != b.nRows() ) {
  249.       error = NEDIM ;
  250.       return 0 ;
  251.    } // if
  252.  
  253.    INDEX  i, nonzero=0, iperm, j ;
  254.    DOUBLE sum, zero = 0.0;
  255.  
  256.    /*************************
  257.     Check for zero diagonals
  258.    **************************/
  259.    for ( i = 1; i <= n ; i++ )
  260.       if ( fabs( A(i,i) ) < tol )
  261.          return i ;
  262.  
  263.    /********************
  264.     forward substitution
  265.    *********************/
  266.  
  267.    for ( i = 1; i <= n ; i++ ) {
  268.       sum = (DOUBLE) b(i) ;
  269.       for ( j = 1; j < i ; j++ )
  270.          sum -= A(j,i) * b(j) ;
  271.       b(i) = (REAL) ( sum / A(i,i) ) ;
  272.    } // for i
  273.  
  274.    /***************************************
  275.     Transform b allowing for trailing zeros
  276.    ****************************************/
  277.  
  278.    for ( i = n ; i >= 1 ; i-- ) {
  279.       sum = (DOUBLE) b(i) ;
  280.       if ( nonzero ) {
  281.          for ( j = i+1 ; j <= nonzero ; j++ )
  282.             sum -= A(j,i) * b(j) ;
  283.       } else if ( sum != zero )
  284.          nonzero = i ;
  285.       iperm = index(i) ;
  286.       b(i) = b( iperm ) ;
  287.       b(iperm) = (REAL) sum ;
  288.    } // for i
  289.    return 0 ;
  290. } // lutSolve
  291.  
  292. INDEX luHager( const matrix& A, const indexArray& index, REAL& est,
  293.            REAL tol, matError& error, INDEX iter )
  294. /****************************************************************
  295.    Estimates lower bound for norm1 of inverse of A.  Assumes A in
  296.    LU form e.g. has been processed by crout. Returns norm estimate
  297.    in est.  iter sets the maximum number of iterations to be used.
  298.    The return value indicates the number of iterations remaining
  299.    on exit from loop, hence if this is non-zero the processed
  300.    "converged".  This routine uses Hager's Convex Optimisation
  301.    Algorithm. See Applied Numerical Linear Algebra, p139 &
  302.    SIAM J Sci Stat Comp 1984 pp 311-16
  303. ****************************************************************/
  304. {
  305.    static char *mName = "luHager" ;
  306.    matFunc func( mName ) ; A.debugInfo( func ) ;
  307.    INDEX i , n, imax ;
  308.    DOUBLE maxz, absz, product, ynorm1, inv_norm1 = 0.0 ;
  309.    INDEX stop ;
  310.    n = A.nRows() ;
  311.    matrix b(n), y(n), z(n) ;
  312.    error = NOERROR ;
  313.    b = (REAL) ( 1.0 / n ) ;
  314.    est = -1.0 ;
  315.    do {
  316.       y = b ;
  317.       if ( luSolve( A, index, y, tol, error ) || error )
  318.          return iter ;
  319.       ynorm1 = y.norm1() ;
  320.       if ( ynorm1 <= inv_norm1 ) {
  321.          stop = TRUE ;
  322.       } else {
  323.          inv_norm1 = ynorm1 ;
  324.          for ( i = 1 ; i <= n ; i++ )
  325.             z(i) = ( y(i) >= 0.0 ? 1.0 : -1.0 ) ;
  326.          if ( lutSolve( A, index, z, tol, error ) || error )
  327.             return iter ;
  328.          imax = 1 ;
  329.          maxz = fabs( (double) z(1) ) ;
  330.          for ( i = 2 ; i <= n ; i++ ) {
  331.             absz = fabs( (double) z(i) ) ;
  332.             if ( absz > maxz ) {
  333.                maxz = absz ;
  334.                imax = i ;
  335.             } // if
  336.          } // for i
  337.          product = (DOUBLE) b.inner(z) ;
  338.          stop = ( maxz <= product ) ;
  339.          if ( !stop ) {
  340.             b = (REAL) 0.0 ;
  341.             b( imax ) = 1.0 ;
  342.          } // if
  343.       } // else
  344.       iter-- ;
  345.    } while ( !stop && iter ) ;
  346.    est = (REAL) inv_norm1 ;
  347.    return iter ;
  348. } // luHager
  349.