home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATDEC.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  11.3 KB  |  430 lines
  1. /**************************************************/
  2. /*      matdec.c source for matDec class          */
  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 "matdec.hpp"
  17.  
  18. /***************************************************/
  19. /*                  matDec methods                */
  20. /***************************************************/
  21.  
  22.  
  23. void matDec::initialValues( void )
  24. {
  25.    status    = 0 ;
  26.    error     = NOERROR ;
  27.    tol       = matrixTol() ;
  28.    det1      = 0 ;
  29.    det2      = 0 ;
  30.    condition = 0 ;
  31. } // matDec initialValues
  32.  
  33. void matDec::copyValues( const matDec& md )
  34. {
  35.    status    = md.status ;
  36.    error     = md.error ;
  37.    tol       = md.tol ;
  38.    det1      = md.det1 ;
  39.    det2      = md.det2 ;
  40.    condition = md.condition ;
  41. } // matDec copyValues
  42.  
  43. #ifdef __cplusplus
  44. matDec::matDec( void ) : matObject()
  45. #else
  46. matDec::matDec( void ) : ()
  47. #endif
  48. {
  49.    initialValues() ;
  50.    m.name( "decM" ) ;
  51. } // matDec( void )
  52.  
  53. #ifdef __cplusplus
  54. matDec::matDec( const matrix& x ) : matObject()
  55. #else
  56. matDec::matDec( const matrix& x ) : ()
  57. #endif
  58. {
  59.    initialValues() ;
  60.    m = x ;
  61.    m.name( "decM" ) ;
  62.    status = ASSIGNED ;
  63.    condition = m.norm1() ;
  64. } // matDec( matrix )
  65.  
  66.  
  67. #ifdef __cplusplus
  68. matDec::matDec( const matDec& md ) : matObject(md)
  69. #else
  70. matDec::matDec( const matDec& md ) : (md)
  71. #endif
  72. {
  73.    m.refer( md.m ) ;
  74.    m.name( "decM" ) ;
  75.    copyValues( md ) ;
  76. } // matDec( matDec& )
  77.  
  78. void matDec::decompose( void )
  79. {
  80.    errorExit( "matDec::decomp", NIMPL ) ;
  81.    return ;
  82. } // matDec::decompose
  83.  
  84. void matDec::solve( matrix& b )
  85. {
  86.    b.errorExit( "matDec::solve", NIMPL ) ;
  87.    return ;
  88. } // matDec::solve
  89.  
  90. void matDec::transSolve( matrix& b )
  91. {
  92.    b.errorExit( "matDec::transSolve", NIMPL ) ;
  93.    return ;
  94. } // matDec::transSolve
  95.  
  96. outFile& matDec::info( outFile& f ) M_CONST
  97. {
  98.    errorExit( "matDec::info", NIMPL ) ;
  99.    return f ;
  100. } // matDec::info
  101.  
  102. void matDec::operator = ( const matDec& md )
  103. {
  104.    static char *mName = "dec op =" ;
  105.    matFunc func( mName ) ; debugInfo( func ) ;
  106.    m = md.m ;
  107.    copyValues( md ) ;
  108. } // matDec = matDec
  109.  
  110. void matDec::assign( const matrix& x )
  111. {
  112.    static char *mName = "dec assign" ;
  113.    matFunc func( mName ) ; // debugInfo( func ) ;
  114.    clear() ;
  115.    m = x ;
  116.    status = ASSIGNED ;
  117.    condition = m.norm1() ;
  118. } // matDec::assign
  119.  
  120. void matDec::capture( matrix& x )
  121. {
  122.    static char *mName = "dec capture" ;
  123.    matFunc func( mName ) ; debugInfo( func ) ;
  124.    initialValues() ;
  125.    m.refer(x) ;
  126.    x.clear() ;
  127.    status = ASSIGNED ;
  128.    condition = m.norm1() ;
  129. } // matDec capture
  130.  
  131. matDec::~matDec( void ) {} // ~matDec
  132.  
  133. void matDec::clear( void )
  134. {
  135.    static char *mName = "dec clear" ;
  136.    matFunc func( mName ) ; debugInfo( func ) ;
  137.    m.clear() ;
  138.    initialValues() ;
  139. } // matDec clear
  140.  
  141. double dfabs( double r )
  142. {
  143.    return ( r < 0.0 ? -r : r ) ;
  144. } // dfabs
  145.  
  146. matError matDec::errorNo( void )
  147. {
  148.    decompose() ;
  149.    return error ;
  150. } // matDec::errorNo
  151.  
  152. INDEX matDec::ok( char *mName )
  153. {
  154.    decompose() ;
  155.    if ( error )
  156.       errorExit( mName, error ) ;
  157.    return OK ;
  158. } // matDec::ok
  159.  
  160. INDEX matDec::ok( void )
  161. {
  162.    decompose() ;
  163.    return !error ;
  164. } // matDec::ok
  165.  
  166. charArray matDec::name( char *newName )
  167. {
  168.    static char *mName = "name" ;
  169.    matFunc func( mName ) ; debugInfo( func ) ;
  170.    if ( newName != 0 )
  171.       nm = newName ;
  172.    return nm ;
  173. } // matDec::name
  174.  
  175. char* matDec::nameStr( void ) M_CONST
  176. {
  177.    return nm.array() ;
  178. } // matDec::name
  179.  
  180. INDEX matDec::Hager( REAL& est, INDEX iter )
  181. /****************************************************************
  182.    Estimates lower bound for norm1 of inverse of A. Returns norm
  183.    estimate in est.  iter sets the maximum number of iterations
  184.    to be used.
  185.    The return value indicates the number of iterations remaining
  186.    on exit from loop, hence if this is non-zero the processed
  187.    "converged".  This routine uses Hager's Convex Optimisation
  188.    Algorithm. See Applied Numerical Linear Algebra, p139 &
  189.    SIAM J Sci Stat Comp 1984 pp 311-16
  190. ****************************************************************/
  191. {
  192.    static char *mName = "luHager" ;
  193.    matFunc func( mName ) ; debugInfo( func ) ;
  194.    INDEX i , n, imax ;
  195.    DOUBLE maxz, absz, product, ynorm1, inv_norm1 = 0.0 ;
  196.    INDEX stop ;
  197.  
  198.    decompose() ;
  199.    if ( error )
  200.       errorExit( mName, error ) ;
  201.    n = m.nRows() ;
  202.    matrix b(n), y(n), z(n) ;
  203.    b = (REAL) ( 1.0 / n ) ;
  204.    est = -1.0 ;
  205.    do {
  206.       y = b ;
  207.       solve(y) ;
  208.       if ( error )
  209.          return iter ;
  210.       ynorm1 = y.norm1() ;
  211.       if ( ynorm1 <= inv_norm1 ) {
  212.          stop = TRUE ;
  213.       } else {
  214.          inv_norm1 = ynorm1 ;
  215.          for ( i = 1 ; i <= n ; i++ )
  216.             z(i) = ( y(i) >= 0.0 ? 1.0 : -1.0 ) ;
  217.          transSolve(z) ;
  218.          if ( error )
  219.             return iter ;
  220.          imax = 1 ;
  221.          maxz = dfabs( (double) z(1) ) ;
  222.          for ( i = 2 ; i <= n ; i++ ) {
  223.             absz = dfabs( (double) z(i) ) ;
  224.             if ( absz > maxz ) {
  225.                maxz = absz ;
  226.                imax = i ;
  227.             } // if
  228.          } // for i
  229.          product = (DOUBLE) b.inner(z) ;
  230.          stop = ( maxz <= product ) ;
  231.          if ( !stop ) {
  232.             b = (REAL) 0.0 ;
  233.             b( imax ) = 1.0 ;
  234.          } // if
  235.       } // else
  236.       iter-- ;
  237.    } while ( !stop && iter ) ;
  238.    est = (REAL) inv_norm1 ;
  239.    return iter ;
  240. } // matDec::Hager
  241.  
  242. REAL matDec::cond( void )
  243. {
  244.    static char *mName = "cond" ;
  245.    matFunc func( mName ) ; debugInfo( func ) ;
  246.    if ( !( status & CONDITION ) ) {
  247.       if ( ok( mName ) ) {
  248.           REAL inorm ;
  249.           if ( Hager( inorm ) && !error )
  250.              condition *= inorm ;
  251.           else if ( error )
  252.              errorExit( mName, error ) ;
  253.           else // no convergence in Hager
  254.              errorExit( mName, NCONV ) ;
  255.       } // if
  256.       status |= CONDITION ;
  257.    } // if
  258.    return condition ;
  259. } // matDec::cond
  260.  
  261. REAL matDec::setTol( const REAL newTol )
  262. {
  263.    REAL oldTol = tol ;
  264.    if ( newTol >= 0.0 )
  265.       tol = newTol ;
  266.    return oldTol ;
  267. } // matDec::setTol
  268.  
  269. outFile& matDec::decInfo( outFile& f, const char* decName ) M_CONST
  270. {
  271.    if ( matListCtrl() > 4 )
  272.       return f ;
  273.    objectInfo( f ) ;
  274.    putName( nameStr(), f ) ;
  275.    putName( decName, f ) ;
  276.    if ( status & ASSIGNED )
  277.      putField( m.identity(), f ) ;
  278.    else
  279.      f.write( "NA" ) ;
  280.    f.newLine() ;
  281.    return f ;
  282. } // matDec::decInfo
  283.  
  284. outFile& matDec::put( outFile& f ) M_CONST
  285. {
  286.    static char *mName = "put" ;
  287.    matFunc func( mName ) ; debugInfo( func ) ;
  288.    f.putIndex( status ).putIndex( (INDEX)(error) ) ;
  289.    f.putReal( tol ).newLine() ;
  290.    f.putReal( det1 ).putReal( det2 ) ;
  291.    f.putReal( condition ).newLine() ;
  292.    f.putIndex( m.nRows() ).putIndex( m.nCols() ).newLine() ;
  293.    m.put(f) ;
  294.    return f ;
  295. } // matDec::put
  296.  
  297. outFile& matDec::print( char* decName, outFile& f ) M_CONST
  298. {
  299.    static char *mName = "put" ;
  300.    matFunc func( mName ) ; debugInfo( func ) ;
  301.    f.write( "Status    : " ).writeIndex( status ).newLine() ;
  302.    f.write( "Error     : " ).writeIndex( (INDEX)(error) ).newLine() ;
  303.    f.write( "Tolerance : " ).writeReal( tol ).newLine() ;
  304.    if ( status & DETERMINED ) {
  305.       f.write( "Det       : " ).writeReal( det1 ) ;
  306.       f.write( "  2^ " ).writeReal( det2 ).newLine() ;
  307.    } // if
  308.    if ( status & CONDITION )
  309.       f.write( "Condition : " ).writeReal( condition ).newLine() ;
  310.    f.newLine() ;
  311.    f.write( decName ).newLine() ;
  312.    INDEX oldFormat = matFormat( DISPLAY ) ;
  313.    m.put(f) ;
  314.    matFormat( oldFormat ) ;
  315.    return f ;
  316. } // matDec::print
  317.  
  318. inFile& matDec::get( inFile& f )
  319. {
  320.    static char *mName = "get" ;
  321.    matFunc func( mName ) ; debugInfo( func ) ;
  322.    INDEX err ;
  323.    f.getIndex( status ).getIndex( err ) ;
  324.    f.getReal( tol ).nextLine() ;
  325.    f.getReal( det1 ).getReal( det2 ) ;
  326.    f.getReal( condition ).nextLine() ;
  327.    error = matError( err ) ;
  328.    INDEX nr, nc ;
  329.    f.getIndex( nr ).getIndex( nc ).nextLine() ;
  330.    m.reset( nr, nc ) ;
  331.    m.get(f) ;
  332.    return f ;
  333. } // matDec::get
  334.  
  335. void matDec::multiSolve( matrix& B )
  336. //   Solve set of equations with RHS in columns of B
  337. {
  338.    static char *mName = "multiSolve" ;
  339.    matFunc func( mName ) ; debugInfo( func ) ;
  340.    if ( m.nRows() != B.nRows() )
  341.       errorExit( mName, NEDIM ) ;
  342.    INDEX n = B.nCols(), i ;
  343.    refMatrix b( B ) ;
  344.    for ( i = 1 ; i <= n && !error ; i++ ) {
  345.       b.refCol(i) ;
  346.       solve( b ) ;
  347.    } // for i
  348.    if ( error )
  349.       errorExit( mName, error ) ;
  350. } // matDec::multiSolve
  351.  
  352. void matDec::inverse( matrix& inv )
  353. {
  354.    static char *mName = "multiSolve" ;
  355.    matFunc func( mName ) ; debugInfo( func ) ;
  356.    INDEX nr = m.nRows(), nc = m.nCols() ;
  357.    if ( inv.nRows() != nr || inv.nCols() != nc )
  358.       inv.reset( nr, nc ) ;
  359.    inv = 0 ;
  360.    inv.setDiag(1.0) ;
  361.    multiSolve( inv ) ;
  362. } // matDec::inverse
  363.  
  364. void dProduct( const matrix& A, REAL& d1, REAL& d2,
  365.            REAL tol, matError& error )
  366. /*********************************************************
  367.    Returns product of A's elements in d1 and d2.
  368.    d1 is a mantissa and d2 an exponent for powers of 2.
  369.    If A is in diagonal of triangular matrix form this will
  370.    be determinant. Assumes A is column vector.
  371.    Based on Bowler, Martin, Peters and Wilkinson in HACLA
  372. ***********************************************************/
  373. {
  374.    REAL  t1   = 1.0, t2  = 0.0 ;
  375.    REAL  zero = 0.0, one = 1.0, four = 4.0, sixteen = 16.0 ;
  376.    REAL  sixteenth = 0.0625 ;
  377.    INDEX n = A.nRows() ;
  378.  
  379.    error = NOERROR ;
  380.    for ( INDEX i = 1 ; ( i <= n ) && ( t1 != zero ) ; i++ ) {
  381.       if ( dfabs( A(i) ) > tol ) {
  382.          t1 *= (DOUBLE) A(i) ;
  383.          while ( dfabs( t1 ) > one ) {
  384.             t1 *= sixteenth ;
  385.             t2 += four ;
  386.          } // while
  387.          while ( dfabs( t1 ) < sixteenth ) {
  388.             t1 *= sixteen ;
  389.             t2 -= four ;
  390.          } // while
  391.       } else {
  392.         t1 = zero ;
  393.         t2 = zero ;
  394.       } // else
  395.    } // for
  396.    d1 = (REAL) t1 ;
  397.    d2 = (REAL) t2 ;
  398.    return ;
  399. } // dProduct
  400.  
  401. void matDec::det( REAL& d1, REAL& d2 )
  402. {
  403.    static char *mName = "det" ;
  404.    matFunc func( mName ) ; debugInfo( func ) ;
  405.    if ( !( status & DETERMINED ) ) {
  406.       decompose() ;
  407.       if ( error == SINGM ) {
  408.          det1 = 0.0 ;
  409.          det2 = 0.0 ;
  410.       } else if ( ok( mName ) ) {
  411.          dProduct( m.diag(), det1, det2, tol, error ) ;
  412.          if ( error )
  413.             errorExit( mName, error ) ;
  414.       } // else
  415.       status |= DETERMINED ;
  416.    } // if
  417.    d1 = det1 ;
  418.    d2 = det2 ;
  419. } // matDec::det
  420.  
  421. void matDec::reportDet( outFile& fout )
  422. {
  423.    static char *mName = "reportDet" ;
  424.    matFunc func( mName ) ; debugInfo( func ) ;
  425.    REAL d1, d2 ;
  426.    det( d1, d2 ) ;
  427.    fout.write( "det = " ).writeReal( d1 ) ;
  428.    fout.write( " x  2^  " ).writeReal( d2 ).newLine( 2 ) ;
  429. } // matDec::reportDet
  430.