home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATSUB.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  6.8 KB  |  225 lines
  1. /**************************************************/
  2. /*     matsub.c source for submatrix 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 "matrix.hpp"
  17.  
  18.  
  19. /*********************************************************/
  20. /*   reference subMatrix methods and functions           */
  21. /*********************************************************/
  22.  
  23.  
  24. void matrix::setSub( INDEX r1, INDEX r2, INDEX c1, INDEX c2,
  25.                      const matrix& y )
  26. {
  27.    static char *mName = "setSub" ;
  28.    matFunc func( mName ) ; debugInfo( func ) ;
  29.    if ( r1 < 1 || r1 > r2 || r2 > nRows() ||
  30.         c1 < 1 || c1 > c2 || c2 > nCols() )
  31.       errorExit( mName, NRANG ) ;
  32.    INDEX nr = r2 - r1 + 1 , nc = c2 - c1 + 1 ;
  33.    if ( nr != y.nRows() || nc != y.nCols() )
  34.       errorExit( mName, NEDIM ) ;
  35.    INDEX i, j, p, q ;
  36.    for ( j = 1, q = c1 ; j <= nc ; j++, q++ )
  37.       for ( i = 1, p = r1 ; i <= nr ; i++, p++ )
  38.      mat(p,q) = y(i,j) ; 
  39. } // matrix::setSub
  40.  
  41. void matrix::setRow( INDEX i, const matrix& y )
  42. {
  43.    static char *mName = "setRow" ;
  44.    matFunc func( mName ) ; debugInfo( func ) ;
  45.    setSub( i, i, 1, nCols(), y ) ;
  46. } // matrix::setRow
  47.  
  48. void matrix::setCol( INDEX j, const matrix& y )
  49. {
  50.    static char *mName = "setCol" ;
  51.    matFunc func( mName ) ; debugInfo( func ) ;
  52.    setSub( 1, nRows(), j, j, y ) ;
  53. } // matrix::setCol
  54.  
  55. void matrix::refer( const matrix& x )
  56. {
  57.    static char *mName = "refer" ;
  58.    matFunc func( mName ) ; debugInfo( func ) ;
  59.    if ( pm != 0 && --pm->ref == 0 )
  60.       delete pm ;
  61.    pm = x.pm ;
  62.    pm->ref++ ;
  63. } // matrix:refer
  64.  
  65. matrix matrix::sub( INDEX r1, INDEX r2, 
  66.                        INDEX c1, INDEX c2  ) M_CONST
  67. {
  68.    static char *mName = "sub" ;
  69.    matFunc func( mName ) ; debugInfo( func ) ;
  70.    matrix z ; 
  71.    z.pm = new matMap( *pm, r1, r2, c1, c2 ) ;
  72.    return z ;
  73. } // matrix::sub
  74.  
  75. matrix matrix::row( INDEX r1, INDEX r2 ) M_CONST
  76. {
  77.    static char *mName = "row" ;
  78.    matFunc func( mName ) ; debugInfo( func ) ;
  79.    if ( r2 == 0 )
  80.       r2 = r1 ;
  81.    return sub( r1, r2, 1, nCols() ) ;
  82. } // matrix::row
  83.  
  84. matrix matrix::col( INDEX c1, INDEX c2 ) M_CONST
  85. {
  86.    static char *mName = "col" ;
  87.    matFunc func( mName ) ; debugInfo( func ) ;
  88.    if ( c2 == 0 )
  89.       c2 = c1 ;
  90.    return sub( 1, nRows(), c1, c2 ) ;
  91. } // matrix::col
  92.  
  93.  
  94. /*********************************************************/
  95. /*        non-reference methods and functions            */
  96. /*********************************************************/
  97.  
  98.  
  99. void matrix::setDiag( matrix& dg, int k )
  100. {
  101.    static char *mName = "setDiag" ;
  102.    matFunc func( mName ) ; debugInfo( func ) ;
  103.    int nr = nRows(), nc = nCols(), nd, offc, offr ;
  104.    if ( dg.nCols() != 1 )
  105.       errorExit( mName, NOTVC ) ;
  106.    if ( ( k > nc - 1 ) || ( nr + k - 1 < 0 ) )
  107.       errorExit( mName, NRANG ) ;
  108.    offr = k < 0 ? -k : 0 ;
  109.    offc = k > 0 ?  k : 0 ;
  110.    nd = (INDEX) ( nr + k <= nc ? nr - offr : nc - offc ) ;
  111.    if ( nd != dg.nRows() )
  112.       errorExit( mName, NEDIM ) ;
  113.    for ( int j = 1 ; j <= nd ; j++ )
  114.       mat( offr + j, offc + j ) = dg(j) ;
  115. } // matrix::setDiag()
  116.  
  117. void matrix::setDiag( REAL r, int k )
  118. {
  119.    static char *mName = "setDiag(r)" ;
  120.    matFunc func( mName ) ; debugInfo( func ) ;
  121.    int nr = nRows(), nc = nCols(), nd, offc, offr ;
  122.    if ( ( k > nc - 1 ) || ( nr + k - 1 < 0 ) )
  123.       errorExit( mName, NRANG ) ;
  124.    offr = k < 0 ? -k : 0 ;
  125.    offc = k > 0 ?  k : 0 ;
  126.    nd = (INDEX) ( nr + k <= nc ? nr - offr : nc - offc ) ;
  127.    for ( int j = 1 ; j <= nd ; j++ )
  128.       mat( offr + j, offc + j ) = r ;
  129. } // matrix::setDiag( REAL )
  130.  
  131. matrix& matrix::diagOf( const matrix& x, int k )
  132. {
  133.    static char *mName = "diagOf" ;
  134.    matFunc func( mName ) ; debugInfo( func ) ;
  135.    int nr = x.nRows(), nc = x.nCols(), nd, offr, offc ;
  136.    if ( ( k > nc - 1 ) || ( nr + k - 1 < 0 ) )
  137.       errorExit( mName, NRANG ) ;
  138.    offr = k < 0 ? -k : 0 ;
  139.    offc = k > 0 ?  k : 0 ;
  140.    nd = (INDEX) ( nr + k <= nc ? nr - offr : nc - offc ) ;
  141.    reset( nd ) ;
  142.    for ( int j = 1 ; j <= nd ; j++ )
  143.       mat(j) = x( offr + j, offc + j ) ;
  144.    return *this ;
  145. } // matrix::diagOf
  146.  
  147. matrix matrix::diag( int k ) M_CONST
  148. {
  149.    static char *mName = "diag" ;
  150.    matFunc func( mName ) ; debugInfo( func ) ;
  151.    matrix z ;
  152.    z.diagOf( *this, k ) ;
  153.    return z ;
  154. } // matrix::diag
  155.  
  156. matrix& matrix::subOf( const matrix& x, INDEX r1, INDEX r2, 
  157.                          INDEX c1, INDEX c2 ) 
  158. {
  159.    static char *mName = "subOf" ;
  160.    matFunc func( mName ) ; debugInfo( func ) ;
  161.    if ( r1 < 1 || r1 > r2 || r2 > x.nRows() || 
  162.         c1 < 1 || c1 > c2 || c2 > x.nCols() )
  163.       errorExit( mName, NRANG ) ;
  164.    INDEX nr = r2 - r1 + 1, nc = c2 - c1 + 1, i, j, p, q ;
  165.    reset( nr, nc ) ;
  166.    for ( j = 1, q = c1 ; j <= nc ; j++, q++ )
  167.       for ( i = 1, p = r1 ; i <= nr ; i++, p++ )
  168.          mat(i,j) = x( p, q ) ;
  169.    return *this ;
  170. } // matrix::subOf
  171.  
  172. matrix& matrix::rowOf( const matrix& x, INDEX i, INDEX j )
  173. {
  174.    static char *mName = "rowOf" ;
  175.    matFunc func( mName ) ; debugInfo( func ) ;
  176.    if ( j == 0 )
  177.       j = i ;
  178.    return subOf( x, i, j, 1, x.nCols() ) ;
  179. } // matrix::rowOf
  180.  
  181. matrix& matrix::colOf( const matrix& x, INDEX i, INDEX j )
  182. {
  183.    static char *mName = "colOf" ;
  184.    matFunc func( mName ) ; debugInfo( func ) ;
  185.    if ( j == 0 )
  186.       j = i ;
  187.    return subOf( x, 1, x.nRows(), i, j ) ;
  188. } // matrix::colOf
  189.  
  190. matrix& matrix::triuOf( const matrix& x, int k )
  191. {
  192.    static char *mName = "triuOf" ;
  193.    matFunc func( mName ) ; debugInfo( func ) ;
  194.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  195.    reset( nr, nc ) ;
  196.    int h ;
  197.    for ( i = 1 ; i <= nr ; i++ ) {
  198.       h = ( i + k <= nc ) ? i + k : nc + 1 ;
  199.       h = ( h > 1 ) ? h : 1 ;
  200.       for ( j = 1 ; j < (INDEX) h ; j++ )
  201.          mat(i,j) = 0.0 ;
  202.       for ( j = (INDEX) h ; j <= nc ; j++ )
  203.          mat(i,j) = x(i,j) ;
  204.    } // for i
  205.    return *this ;
  206. } // matrix triuOf
  207.  
  208. matrix& matrix::trilOf( const matrix& x, int k )
  209. {
  210.    static char *mName = "trilOf" ;
  211.    matFunc func( mName ) ; debugInfo( func ) ;
  212.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  213.    reset( nr, nc ) ;
  214.    int h ;
  215.    for ( i = 1 ; i <= nr ; i++ ) {
  216.       h = ( i + k < 1 ) ? 0 : i + k ;
  217.       h = ( h <= nc ) ? h : nc ;
  218.       for ( j = 1 ; j <= (INDEX) h ; j++ )
  219.          mat(i,j) = x(i,j) ;
  220.       for ( j = (INDEX) h + 1 ; j <= nc ; j++ )
  221.          mat(i,j) = 0.0 ;
  222.    } // for i
  223.    return *this ;
  224. } // matrix trilOf
  225.