home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATGEN.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  19.1 KB  |  676 lines
  1. /**************************************************/
  2. /* matgen.c source for general MatClass 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. #include <math.h>
  19.  
  20. matrix& matrix::linear( REAL a, const matrix& x,
  21.                         const matrix& y )
  22. /***********************************************
  23.          return   a * x + y  in this
  24. ************************************************/
  25. {
  26.    static char *mName = "linear(a,x,y)" ;
  27.    matFunc func( mName ) ; debugInfo( func ) ;
  28.    x.checkDims( y, mName ) ;
  29.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  30.    reset( nr, nc ) ;
  31.    if ( a == 1.0 ) {
  32.       for ( j = 1; j <= nc ; j++ )
  33.          for ( i = 1; i <= nr ; i++ )
  34.             mat(i,j) = x(i,j) + y(i,j) ;
  35.    } else if ( a == -1.0 ) {
  36.       for ( j = 1; j <= nc ; j++ )
  37.          for ( i = 1; i <= nr ; i++ )
  38.             mat(i,j) = y(i,j) - x(i,j) ;
  39.    } else if ( a == 0.0 ) {
  40.       for ( j = 1; j <= nc ; j++ )
  41.          for ( i = 1; i <= nr ; i++ )
  42.             mat(i,j) = y(i,j) ;
  43.    } else {
  44.       for ( j = 1; j <= nc ; j++ )
  45.          for ( i = 1; i <= nr ; i++ )
  46.             mat(i,j) = a * x(i,j) + y(i,j) ;
  47.    } // else
  48.    return *this ;
  49. } // matrix::linear(a,x,y)
  50.  
  51. matrix& matrix::linear( REAL a, const matrix& x,
  52.                         REAL b )
  53. /***********************************************
  54.         return  a * x + b  in  this
  55. ***********************************************/
  56. {
  57.    static char *mName = "linear(a,x,b)" ;
  58.    matFunc func( mName ) ; debugInfo( func ) ;
  59.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  60.    reset( nr, nc ) ;
  61.    if ( a == 1.0 ) {
  62.       for ( j = 1; j <= nc ; j++ )
  63.          for ( i = 1; i <= nr ; i++ )
  64.             mat(i,j) = x(i,j) + b ;
  65.    } else if ( a == -1.0 ) {
  66.       for ( j = 1; j <= nc ; j++ )
  67.          for ( i = 1; i <= nr ; i++ )
  68.             mat(i,j) = b - x(i,j) ;
  69.    } else if ( a == 0.0 ) {
  70.       for ( j = 1; j <= nc ; j++ )
  71.          for ( i = 1; i <= nr ; i++ )
  72.             mat(i,j) = b ;
  73.    } else {
  74.       for ( j = 1; j <= nc ; j++ )
  75.          for ( i = 1; i <= nr ; i++ )
  76.             mat(i,j) = a * x(i,j) + b ;
  77.    } // else
  78.    return *this ;
  79. } // linear(a,x,b)
  80.  
  81. matrix& matrix::minus( void )
  82. {
  83.    static char *mName = "minus" ;
  84.    matFunc func( mName ) ; debugInfo( func ) ;
  85.    INDEX nr = nRows(), nc = nCols(), i, j ;
  86.    for ( j = 1; j <= nc ; j++ )
  87.       for ( i = 1; i <= nr ; i++ )
  88.          mat(i,j) = - mat(i,j) ;
  89.    return *this ;
  90. } // matrix minus
  91.  
  92. matrix& matrix::operator += ( const matrix &y )
  93. {
  94.    static char *mName = "op +=" ;
  95.    matFunc func( mName ) ; debugInfo( func ) ;
  96.    INDEX nr = nRows(), nc = nCols(), i, j ;
  97.    for ( j = 1; j <= nc ; j++ )
  98.       for ( i = 1; i <= nr ; i++ )
  99.          mat(i,j) += y(i,j) ;
  100.    return *this ;
  101. } // matrix +=
  102.  
  103. matrix& matrix::operator -= ( const matrix& y )
  104. {
  105.    static char *mName = "op -=" ;
  106.    matFunc func( mName ) ; debugInfo( func ) ;
  107.    INDEX nr = nRows(), nc = nCols(), i, j ;
  108.    for ( j = 1; j <= nc ; j++ )
  109.       for ( i = 1; i <= nr ; i++ )
  110.          mat(i,j) -= y(i,j) ;
  111.    return *this ;
  112. } // matrix op -=
  113.  
  114. matrix operator - ( const matrix& x )
  115. {
  116.    static char *mName = "unary -" ;
  117.    matFunc func( mName ) ; x.debugInfo( func ) ;
  118.    matrix z ;
  119.    z.linear( -1.0, x, 0.0 ) ;
  120.    return z ;
  121. } // matrix unary -
  122.  
  123. matrix operator + ( const matrix& x, const matrix &y )
  124. {
  125.    static char *mName = "op +" ;
  126.    matFunc func( mName ) ; x.debugInfo( func ) ;
  127.    matrix z ;
  128.    z.linear( 1.0, x, y ) ;
  129.    return z ;
  130. } // matrix op +
  131.  
  132. matrix operator - ( const matrix& x, const matrix &y )
  133. {
  134.    static char *mName = "op -" ;
  135.    matFunc func( mName ) ; x.debugInfo( func ) ;
  136.    matrix z ;
  137.    z.linear( -1.0, y, x ) ;
  138.    return z ;
  139. } // matrix op -
  140.  
  141.  
  142. matrix& matrix::operator  = ( REAL r )
  143. {
  144.    static char *mName = "matrix=REAL" ;
  145.    matFunc func( mName ) ; debugInfo( func ) ;
  146.    INDEX nr = nRows(), nc = nCols(), i, j ;
  147.    for ( j = 1 ; j <= nc ; j++ )
  148.       for ( i = 1 ; i <= nr ; i++ )
  149.          mat(i,j) = r ;
  150.    return *this ;
  151. } // matrix = REAL
  152.  
  153. matrix& matrix:: operator += ( REAL r )
  154. {
  155.    static char *mName = "matrix+=REAL" ;
  156.    matFunc func( mName ) ; debugInfo( func ) ;
  157.    INDEX nr = nRows(), nc = nCols(), i, j ;
  158.    for ( j = 1 ; j <= nc ; j++ )
  159.       for ( i = 1 ; i <= nr ; i++ )
  160.          mat(i,j) += r ;
  161.    return *this ;
  162. } // matrix += REAL
  163.  
  164. matrix& matrix::operator -= ( REAL r )
  165. {
  166.    static char *mName = "matrix-=REAL" ;
  167.    matFunc func( mName ) ; debugInfo( func ) ;
  168.    *this += (-r) ;
  169.    return *this ;
  170. } // matrix -= REAL
  171.  
  172. matrix& matrix::operator *= ( REAL r )
  173. {
  174.    static char *mName = "matrix*=REAL" ;
  175.    matFunc func( mName ) ; debugInfo( func ) ;
  176.    INDEX nr = nRows(), nc = nCols(), i, j ;
  177.    for ( j = 1 ; j <= nc ; j++ )
  178.       for ( i = 1 ; i <= nr ; i++ )
  179.          mat(i,j) *= r ;
  180.    return *this ;
  181. } // matrix *= REAL
  182.  
  183. matrix& matrix::operator /= ( REAL r )
  184. {
  185.    static char *mName = "matrix/=REAL" ;
  186.    matFunc func( mName ) ; debugInfo( func ) ;
  187.    if ( r == 0.0 )
  188.       errorExit( mName, ZEROD ) ;
  189.    *this *= (1/r) ;
  190.    return *this ;
  191. } // matrix *= REAL
  192.  
  193. matrix operator + ( const matrix& x, REAL r )
  194. {
  195.    static char *mName = "matrix+REAL" ;
  196.    matFunc func( mName ) ; x.debugInfo( func ) ;
  197.    matrix z ;
  198.    z.linear( 1.0, x, r ) ;
  199.    return z ;
  200. } // matrix + REAL
  201.  
  202. matrix operator * ( const matrix& x, REAL r )
  203. {
  204.    static char *mName = "matrix*REAL" ;
  205.    matFunc func( mName ) ; x.debugInfo( func ) ;
  206.    matrix z ;
  207.    z.linear( r, x, 0.0 ) ;
  208.    return z ;
  209. } // matrix * REAL
  210.  
  211. matrix operator - ( const matrix& x, REAL r )
  212. {
  213.    static char *mName = "matrix-REAL" ;
  214.    matFunc func( mName ) ; x.debugInfo( func ) ;
  215.    matrix z ;
  216.    z.linear( 1.0, x, -r ) ;
  217.    return z ;
  218. } // matrix - REAL
  219.  
  220. matrix operator / ( const matrix& x, REAL r )
  221. {
  222.    static char *mName = "matrix/REAL" ;
  223.    matFunc func( mName ) ; x.debugInfo( func ) ;
  224.    if ( r == 0.0 )
  225.       x.errorExit( mName, ZEROD ) ;
  226.    matrix z ;
  227.    z.linear( 1.0 / r, x, 0.0 ) ;
  228.    return z ;
  229. } // matrix - REAL
  230.  
  231. matrix operator + ( REAL r, const matrix &y )
  232. {
  233.    static char *mName = "REAL+matrix" ;
  234.    matFunc func( mName ) ; y.debugInfo( func ) ;
  235.    matrix z ;
  236.    z.linear( 1.0, y, r ) ;
  237.    return z ;
  238. } // REAL + matrix
  239.  
  240. matrix operator - ( REAL r, const matrix &y )
  241. {
  242.    static char *mName = "REAL-matrix" ;
  243.    matFunc func( mName ) ; y.debugInfo( func ) ;
  244.    matrix z ;
  245.    z.linear( -1.0, y, r ) ;
  246.    return z ;
  247. } // REAL - matrix
  248.  
  249. matrix operator * ( REAL r, const matrix &y )
  250. {
  251.    static char *mName = "REAL*matrix" ;
  252.    matFunc func( mName ) ; y.debugInfo( func ) ;
  253.    matrix z ;
  254.    z.linear( r, y, 0.0 ) ;
  255.    return z ;
  256. } // REAL * matrix
  257.  
  258. matrix operator / ( REAL r, const matrix& x )
  259. /*
  260.    Element wise division of r by x
  261. */
  262. {
  263.    static char *mName = "REAL/matrix" ;
  264.    matFunc func( mName ) ; x.debugInfo( func ) ;
  265.    REAL s ;
  266.    INDEX nr = x.nRows(), nc = x.nCols() ;
  267.    matrix z( nr, nc ) ;
  268.    if ( x.isNull() )
  269.       x.errorExit( mName, NULRF ) ;
  270.    for ( INDEX j = 1 ; j <= nc ; j++ ) {
  271.       for ( INDEX i = 1 ; i <= nr ; i++ ) {
  272.          if ( ( s = x(i,j) ) == 0.0 ) {
  273.             x.errorij( i, j ) ;
  274.             x.errorExit( mName, ZEROD ) ;
  275.          } // if
  276.          z.mat(i,j) = r / s ;
  277.       } // for j
  278.    } // for i
  279.    return z ;
  280. } // REAL / matrix
  281.  
  282. matrix& matrix::multOf( const matrix& x,
  283.                         const matrix& y )
  284. /******************************************
  285.      Return multiple of x and y in this
  286. ******************************************/
  287. {
  288.    static char* mName = "multOf" ;
  289.    matFunc func( mName ) ; debugInfo( func ) ;
  290.    INDEX n, nr, nc, i, j, k ;
  291.    DOUBLE r ;
  292.    if ( x.nCols() != y.nRows() )
  293.       x.errorExit( mName, NEDIM ) ;
  294.    n  = x.nCols() ;
  295.    nr = x.nRows() ;
  296.    nc = y.nCols() ;
  297.    reset( nr, nc ) ;
  298.    for ( i = 1 ; i <= nr ; i++ ) {
  299.       for ( j = 1 ; j <= nc ; j++ ) {
  300.          r = 0.0 ;
  301.          for ( k = 1 ; k <= n ; k++ )
  302.             r += x(i,k) * y(k,j) ;
  303.          mat(i,j) =  (REAL) r ;
  304.       } // for j
  305.    } // for i
  306.    return *this ;
  307. } // matrix multOf
  308.  
  309. matrix& matrix::TMultOf( const matrix& x,
  310.                          const matrix& y )
  311. /*******************************************
  312.     Return multiple of x' and y in this
  313. *******************************************/
  314. {
  315.    static char* mName = "TMultOf" ;
  316.    matFunc func( mName ) ; debugInfo( func ) ;
  317.    INDEX n, nr, nc, i, j, k ;
  318.    DOUBLE r ;
  319.    if ( x.nRows() != y.nRows() )
  320.       x.errorExit( mName, NEDIM ) ;
  321.    n  = x.nRows() ;
  322.    nr = x.nCols() ;
  323.    nc = y.nCols() ;
  324.    reset( nr, nc ) ;
  325.    for ( i = 1 ; i <= nr ; i++ ) {
  326.       for ( j = 1 ; j <= nc ; j++ ) {
  327.          r = 0.0 ;
  328.          for ( k = 1 ; k <= n ; k++ )
  329.             r += x(k,i) * y(k,j) ;
  330.          mat(i,j) =  (REAL) r ;
  331.       } // for j
  332.    } // for i
  333.    return *this ;
  334. } // TMultOf
  335.  
  336. matrix matrix::TMult( const matrix& y ) M_CONST
  337. {
  338.    static char* mName = "TMult" ;
  339.    matFunc func( mName ) ; debugInfo( func ) ;
  340.    matrix z ;
  341.    z.TMultOf( *this, y ) ;
  342.    return z ;
  343. } // matrix::TMult
  344.  
  345. matrix& matrix::multTOf( const matrix& x,
  346.                          const matrix& y )
  347. /****************************************
  348.    Return multiple of x and y' in z
  349. ****************************************/
  350. {
  351.    static char* mName = "multTOf" ;
  352.    matFunc func( mName ) ; debugInfo( func ) ;
  353.    INDEX n, nr, nc, i, j, k ;
  354.    DOUBLE r ;
  355.    if ( x.nCols() != y.nCols() )
  356.       x.errorExit( mName, NEDIM ) ;
  357.    n  = x.nCols() ;
  358.    nr = x.nRows() ;
  359.    nc = y.nRows() ;
  360.    reset( nr, nc ) ;
  361.    for ( i = 1 ; i <= nr ; i++ ) {
  362.       for ( j = 1 ; j <= nc ; j++ ) {
  363.          r = 0.0 ;
  364.          for ( k = 1 ; k <= n ; k++ )
  365.             r += x(i,k) * y(j,k) ;
  366.          mat(i,j) =  (REAL) r ;
  367.       } // for j
  368.    } // for i
  369.    return *this ;
  370. } // multTOf
  371.  
  372. matrix matrix::multT( const matrix& y ) M_CONST
  373. {
  374.    static char* mName = "multT" ;
  375.    matFunc func( mName ) ; debugInfo( func ) ;
  376.    matrix z ;
  377.    z.multTOf( *this, y ) ;
  378.    return z ;
  379. } // matrix::multT
  380.  
  381. matrix& matrix::TMultTOf( const matrix& x,
  382.                           const matrix& y )
  383. /******************************************
  384.    Return multiple of x' and y' in z
  385. ******************************************/
  386. {
  387.    static char* mName = "TMultTOf" ;
  388.    matFunc func( mName ) ; debugInfo( func ) ;
  389.    INDEX n, nr, nc, i, j, k ;
  390.    DOUBLE r ;
  391.    if ( x.nRows() != y.nCols() )
  392.       x.errorExit( mName, NEDIM ) ;
  393.    n  = x.nRows() ;
  394.    nr = x.nCols() ;
  395.    nc = y.nRows() ;
  396.    reset( nr, nc ) ;
  397.    for ( i = 1 ; i <= nr ; i++ ) {
  398.       for ( j = 1 ; j <= nc ; j++ ) {
  399.          r = 0.0 ;
  400.          for ( k = 1 ; k <= n ; k++ )
  401.             r += x(k,i) * y(j,k) ;
  402.          mat(i,j) =  (REAL) r ;
  403.       } // for j
  404.    } // for i
  405.    return *this ;
  406. } // TMultTOf
  407.  
  408. matrix matrix::TMultT( const matrix& y ) M_CONST
  409. {
  410.    static char* mName = "TMultT" ;
  411.    matFunc func( mName ) ; debugInfo( func ) ;
  412.    matrix z ;
  413.    z.TMultTOf( *this, y ) ;
  414.    return z ;
  415. } // matrix::TMultT
  416.  
  417. matrix operator * ( const matrix& x,
  418.                     const matrix &y )
  419. {
  420.    static char* mName = "op *" ;
  421.    matFunc func( mName ) ; x.debugInfo( func ) ;
  422.    matrix z ;
  423.    if ( x.isTrans() ) {
  424.       if ( y.isTrans() )
  425.          z.TMultTOf(x,y) ;
  426.       else
  427.          z.TMultOf(x,y) ;
  428.    } else if ( y.isTrans() )
  429.       z.multTOf(x,y) ;
  430.    else
  431.       z.multOf( x, y ) ;
  432.    return z ;
  433. } // matrix operator *
  434.  
  435. matrix matrix::multij( const matrix& y ) M_CONST
  436. /***********************************************
  437.    Element wise multiplication of this by y
  438. ***********************************************/
  439. {
  440.    static char *mName = "multijOf" ;
  441.    matFunc func( mName ) ; debugInfo( func ) ;
  442.    checkDims( y, mName  );
  443.    INDEX nr = nRows(), nc = nCols() ;
  444.    matrix z( nr, nc ) ;
  445.    for ( INDEX j = 1 ; j <= nc ; j++ )
  446.       for ( INDEX i = 1 ; i <= nr ; i++ )
  447.          z(i,j) = mat(i,j) * y(i,j) ;
  448.    return z ;
  449. } // matrix::multij
  450.  
  451. matrix& matrix::multijEq( const matrix& y )
  452. {
  453.    static char* mName = "multijEq" ;
  454.    matFunc func( mName ) ; debugInfo( func ) ;
  455.    checkDims( y, mName  );
  456.    INDEX nr = nRows(), nc = nCols() ;
  457.    for ( INDEX j = 1 ; j <= nc ; j++ )
  458.       for ( INDEX i = 1 ; i <= nr ; i++ )
  459.          mat(i,j) *= y(i,j) ;
  460.    return *this ;
  461. } // matrix::multijEq
  462.  
  463. matrix matrix::divij( const matrix& y ) M_CONST
  464. /**********************************************
  465.       Element wise division of this by y
  466. **********************************************/
  467. {
  468.    char *mName = "divijOf" ;
  469.    matFunc func( mName ) ; debugInfo( func ) ;
  470.    checkDims( y, mName  ) ;
  471.    INDEX nr = nRows(), nc = nCols() ;
  472.    matrix z( nr, nc ) ;
  473.    REAL r ;
  474.    for ( INDEX j = 1 ; j <= nc ; j++ ) {
  475.       for ( INDEX i = 1 ; i <= nr ; i++ ) {
  476.          if ( ( r = y(i,j) ) == 0.0 ) {
  477.             y.errorij( i, j ) ;
  478.             y.errorExit( mName, ZEROD ) ;
  479.          } // if
  480.          z(i,j) = mat(i,j) / r ;
  481.       } // for j
  482.    } // for i
  483.    return z ;
  484. } // matrix::divij
  485.  
  486. matrix& matrix::divijEq( const matrix& y )
  487. {
  488.    static char* mName = "divijEq" ;
  489.    matFunc func( mName ) ; debugInfo( func ) ;
  490.    checkDims( y, mName  ) ;
  491.    INDEX nr = nRows(), nc = nCols() ;
  492.    REAL r ;
  493.    for ( INDEX j = 1 ; j <= nc ; j++ ) {
  494.       for ( INDEX i = 1 ; i <= nr ; i++ ) {
  495.          if ( ( r = y(i,j) ) == 0.0 ) {
  496.             y.errorij( i, j ) ;
  497.             y.errorExit( mName, ZEROD ) ;
  498.          } // if
  499.          mat(i,j) /= r ;
  500.       } // for j
  501.    } // for i
  502.    return *this ;
  503. } // matrix::divijEq
  504.  
  505. REAL matrix::inner( const matrix &y ) M_CONST
  506. /********************************************
  507.        Inner product of this and y
  508. ********************************************/
  509. {
  510.    char *mName = "inner" ;
  511.    matFunc func( mName ) ; debugInfo( func ) ;
  512.    checkDims( y, mName ) ;
  513.    INDEX nc = nCols(), nr = nRows() ;
  514.    if ( nc != y.nCols() || nr != y.nRows() )
  515.       errorExit( mName, NEDIM ) ;
  516.    DOUBLE sum = 0.0 ;
  517.    for ( INDEX j = 1; j <= nc ; j++ ) {
  518.       for ( INDEX i = 1; i <= nr ; i++ )
  519.          sum  += mat(i,j) * y(i,j) ;
  520.    } // for j
  521.    return (REAL) sum ;
  522. } // REAL matrix::inner
  523.  
  524. matrix& matrix::rowMultOf( const matrix& x,
  525.                            const matrix& dg )
  526. /*********************************************
  527.   Multiply the rows of x by the elements of
  528.   column dg, this effectively treats dg as
  529.   a diagonal matrix and forms the pre-
  530.   product diag(dg) * x.
  531. *********************************************/
  532. {
  533.    static char* mName = "rowMultOf" ;
  534.    matFunc func( mName ) ; debugInfo( func ) ;
  535.    INDEX nr = nRows(), nc = nCols(), i, j ;
  536.    reset( nr, nc ) ;
  537.    REAL r ;
  538.    if ( dg.nCols() != 1 )
  539.       dg.errorExit( mName, NOTVC ) ;
  540.    if ( dg.nRows() != nr )
  541.       dg.errorExit( mName, NEDIM ) ;
  542.    for ( i = 1 ; i <= nr ; i++ ) {
  543.       r = dg(i) ;
  544.       for ( j = 1 ; j <= nc ; j++ )
  545.          mat(i,j) = r * x(i,j) ;
  546.    } // for i
  547.    return *this ;
  548. } // matrix::rowMultOf
  549.  
  550. matrix matrix::rowMult( const matrix& dg ) M_CONST
  551. /*************************************************
  552.    Multiply the rows of this by the elements of
  553.    column dg, this effectively treats dg as a
  554.    diagonal matrix and forms the pre-product
  555.    diag(dg) * (*this).
  556. *************************************************/
  557. {
  558.    static char* mName = "rowMult" ;
  559.    matFunc func( mName ) ; debugInfo( func ) ;
  560.    INDEX nr = nRows(), nc = nCols() ;
  561.    matrix z( nr, nc ) ;
  562.    z.rowMultOf( *this, dg ) ;
  563.    return z ;
  564. } // matrix::rowMult
  565.  
  566. matrix& matrix::colMultOf( const matrix& x,
  567.                            const matrix &dg )
  568. /***********************************************
  569.   Multiply the columns of this by the elements
  570.   of column dg, this effectively treats dg as
  571.   a diagonal matrix and forms the post-product
  572.   x * diag(dg).
  573. ***********************************************/
  574. {
  575.    static char *mName = "colMultOf" ;
  576.    matFunc func( mName ) ; debugInfo( func ) ;
  577.    INDEX nr = nRows(), nc = nCols(), i, j ;
  578.    reset( nr, nc ) ;
  579.    REAL r ;
  580.    if ( dg.nCols() != 1 )
  581.       dg.errorExit( mName, NOTVC ) ;
  582.    if ( dg.nRows() != nc )
  583.       dg.errorExit( mName, NEDIM ) ;
  584.    for ( j = 1 ; j <= nc ; j++ ){
  585.       r = dg(j) ;
  586.       for ( i = 1 ; i <= nr ; i++ )
  587.          mat(i,j) = r * x(i,j) ;
  588.    } // for j
  589.    return *this ;
  590. } // matrix::colMultOf
  591.  
  592. matrix matrix::colMult( const matrix &dg ) M_CONST
  593. /*************************************************
  594.    Multiply the columns of this by the elements
  595.    of column dg, this effectively treats dg as
  596.    a diagonal matrix and forms the post-product
  597.    (*this) * diag(dg).
  598. *************************************************/
  599. {
  600.    static char *mName = "colMultOf" ;
  601.    matFunc func( mName ) ; debugInfo( func ) ;
  602.    INDEX nr = nRows(), nc = nCols() ;
  603.    matrix z( nr, nc ) ;
  604.    z.colMultOf( *this, dg ) ;
  605.    return z ;
  606. } // matrix::colMult
  607.  
  608. #include <math.h>
  609.  
  610. matrix& matrix::step( REAL first, REAL increment )
  611. {
  612.    char *mName = "step" ;
  613.    matFunc func( mName ) ; debugInfo( func ) ;
  614.    INDEX nr = nRows(), nc = nCols(), i, j ;
  615.    REAL r = first ;
  616.    for ( j = 1 ; j <= nc ; j++ )
  617.       for (  i = 1; i <= nr; i++ ) {
  618.          mat(i,j) = r ;
  619.          r += increment ;
  620.        } // for
  621.    return *this ;
  622. } // step matrix
  623.  
  624. matrix& joinOf( matrix& z, const matrix& x,
  625.                 const matrix& y )
  626. {
  627.    static char *mName = "joinOf" ;
  628.    matFunc func( mName ) ; x.debugInfo( func ) ;
  629.    INDEX xnr = x.nRows(), xnc = x.nCols() ;
  630.    INDEX ynr = y.nRows(), ync = y.nCols(), nc ;
  631.    if ( ynr != xnr )
  632.       x.errorExit( mName, NEDIM ) ;
  633.    nc = xnc + ync ;
  634.    z.reset( xnr, nc ) ;
  635.    z.setSub( 1, xnr, 1, xnc, x ) ;
  636.    z.setSub( 1, xnr, xnc+1, nc, y ) ;
  637.    return z ;
  638. } // joinOf
  639.  
  640. matrix& stackOf( matrix& z, const matrix& x,
  641.                  const matrix& y )
  642. {
  643.    static char *mName = "stackOf" ;
  644.    matFunc func( mName ) ; x.debugInfo( func ) ;
  645.    INDEX xnr = x.nRows(), xnc = x.nCols() ;
  646.    INDEX ynr = y.nRows(), ync = y.nCols(), nr ;
  647.    if ( ync != xnc )
  648.       x.errorExit( mName, NEDIM ) ;
  649.    nr = xnr + ynr ;
  650.    z.reset( nr, xnc ) ;
  651.    z.setSub( 1, xnr, 1, xnc, x ) ;
  652.    z.setSub( xnr+1, nr, 1, xnc, y ) ;
  653.    return z ;
  654. } // stackOf
  655.  
  656. matrix operator | ( const matrix& x, const matrix &y )
  657. {
  658.    static char *mName = "op |" ;
  659.    matFunc func( mName ) ; x.debugInfo( func ) ;
  660.    INDEX nr = x.nRows(), nc = x.nCols() + y.nCols() ;
  661.    matrix z( nr, nc ) ;
  662.    joinOf( z, x, y ) ;
  663.    return z ;
  664. } // matrix op |
  665.  
  666. matrix operator || ( const matrix& x, const matrix &y )
  667. {
  668.    static char *mName = "op ||" ;
  669.    matFunc func( mName ) ; x.debugInfo( func ) ;
  670.    INDEX nr = x.nRows() + y.nRows(), nc = x.nCols() ;
  671.    matrix z( nr, nc ) ;
  672.    stackOf( z, x, y ) ;
  673.    return z ;
  674. } // matrix op ||
  675.  
  676.