home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATQRH.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  12.0 KB  |  368 lines
  1. /**************************************************/
  2. /*        matqrh.c source for matQrh 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 "matqrh.hpp"
  17. #include <math.h>
  18.  
  19. /****************************************************/
  20. /*                 qrhDec methods                   */
  21. /****************************************************/
  22.  
  23. INDEX hReflector( matrix& v, REAL& alpha, REAL& beta,
  24.                   REAL& vMax, REAL tol, matError& error )
  25. /***************************************************************
  26.    Reflector is I - beta * u * u' where
  27.  
  28.                   u = w + alpha * e(1)
  29.                   w = v / max(abs(v)).
  30.               alpha = sign(v(1)) ||w||
  31.                beta = 1 / ( alpha^2 + alpha * w(1) )
  32.                     = 1 / ( alpha * u(1) ).
  33.  
  34.    The use of the scaled w avoids overflow & underflow.
  35.    See Golub and Loan first edition pp 38-41.
  36.    v is overwritten with u. max(abs(v)) is returned in vMax,
  37.    if this is less than tol then proc returns singularity
  38.    error. Proc insists v is a column vector.
  39. ***************************************************************/
  40. {
  41.    static char *mName = "hReflector" ;
  42.    matFunc func( mName ) ; v.debugInfo( func ) ;
  43.    DOUBLE sum ;
  44.    REAL s ;
  45.    error = NOERROR ;
  46.    if ( v.nCols() != 1 ) {
  47.       error = NOTVC ;
  48.       return 0 ;
  49.    } // if
  50.    INDEX i, n = v.nRows() ;
  51.    vMax = REAL( fabs( double( v(1) ) ) ) ;
  52.    for ( i = 2 ; i <= n ; i++ ) {
  53.        s = REAL( fabs( double( v(i) ) ) ) ;
  54.        if ( s > vMax )
  55.           vMax = s;
  56.    } // for i
  57.    if ( vMax < tol ){
  58.       error = SINGM ;
  59.       return 0 ;
  60.    } // if
  61.    // overwrite v with w and form alpha
  62.    sum = 0.0 ;
  63.    for ( i = 1 ; i <= n ; i++ ) {
  64.       v(i) /= vMax ;
  65.       sum += v(i) * v(i) ;
  66.    } // for
  67.    alpha = REAL( sqrt( sum ) ) ;
  68.    if ( v(1) < 0.0 )
  69.       alpha = - alpha ;
  70.    v(1) += alpha ;
  71.    beta = 1.0 / ( alpha * v(1) ) ;
  72.    return OK ;
  73. } // hReflector
  74.  
  75. INDEX hReflect( matrix& y, const matrix& u, REAL beta,
  76.                 matError& error )
  77. /**************************************************
  78.    Reflect y in plane with normal u. Assume
  79.    beta contains 2/u'u. Assumes y and u has same
  80.    number of rows and that u is a column vector.
  81.    u and y must have the same number of rows and
  82.    u must be a column vector.
  83. **************************************************/
  84. {
  85.    static char *mName = "hReflect" ;
  86.    matFunc func( mName ) ; y.debugInfo( func ) ;
  87.    error = NOERROR ;
  88.    if ( u.nCols() != 1 ) {
  89.       error = NOTVC ;
  90.       return 0 ;
  91.    } // if
  92.    if ( u.nRows() != y.nRows() ) {
  93.       error = NEDIM ;
  94.       return 0 ;
  95.    } // if
  96.    INDEX j, n = y.nCols() ;
  97.    refMatrix yCol( y ) ;
  98.    for ( j = 1; j <= n ; j++ ) {
  99.       // yCol reference to jth col of y
  100.       yCol.refCol( j ) ;
  101.       // overwrite column with column - beta * (u'column) * u
  102.       yCol.linear( - beta * u.inner( yCol ), u, yCol ) ;
  103.    } // for j
  104.    return OK ;
  105. } // hReflect
  106.  
  107. INDEX qrh( matrix& x, matrix& diagR, matrix& b, REAL tol,
  108.            matError& error )
  109. /****************************************************************
  110.    QR decomposition of x by Householder transformations. See
  111.    Golub & Loan first edition p41 & Sec 6.2. R is returned in
  112.    upper triang of x and diagR.  Q returned in `u-form' in
  113.    lower triang of x and b, the latter containing the
  114.    "Householder betas". diagR and b are reset if required.
  115.    Errors arise from formation of reflectors i.e. singularity,
  116.    and the reflections, though the latter should not lead to
  117.    difficulties. Note it attempts to handle the cases where
  118.    the number of rows is less than or equal to the number of
  119.    columns.
  120. ****************************************************************/
  121. {
  122.    static char *mName = "qrh" ;
  123.    matFunc func( mName ) ; x.debugInfo( func ) ;
  124.    INDEX nr = x.nRows(), nc = x.nCols(), n, k ;
  125.    refMatrix kCol(x), subMat(x) ;
  126.    REAL colMax, alpha, beta ;
  127.    if ( nr <= nc ) {
  128.       n = nr - 1 ;
  129.       diagR.reset( nr ) ;
  130.       b.reset( nr ) ;
  131.    } else {
  132.       n = nc ;
  133.       diagR.reset( nc ) ;
  134.       b.reset( nc ) ;
  135.    } // else
  136.    for ( k = 1 ; k <= n ; k++ ) {
  137.       kCol.refSub( k, nr, k, k ) ;
  138.       if ( !hReflector( kCol, alpha, beta, colMax, tol, error ) )
  139.          return 0 ;
  140.       diagR(k) = - alpha * colMax ;
  141.       b(k) = beta ;
  142.       if ( k < nc ) {
  143.          // set subMat to sub matrix to right of kCol
  144.          subMat.refSub( k, nr, k+1, nc ) ;
  145.          // apply Householder reflection to subMat
  146.          if ( !hReflect( subMat, kCol, beta, error ) )
  147.             return 0 ;
  148.       } // if
  149.    } // for k
  150.    if ( nr <= nc ) {
  151.       diagR( nr ) = x( nr, nr ) ;
  152.       b( nr ) = 0 ;
  153.    } // if
  154.    return 1 ;
  155. } // qrh
  156.  
  157. INDEX backSub( matrix& y, const matrix& x, REAL tol,
  158.                matError& error )
  159. /*********************************************
  160.    Solve xb = y assuming x is upper triang
  161.    Assume number of equations equals number
  162.    of columns in x. This allows x to have
  163.    more rows than columns e.g. after QR
  164.    decomposition. Proc requires the number
  165.    of rows be at least as great as the
  166.    number of equations.
  167. *********************************************/
  168. {
  169.    static char *mName = "backSub" ;
  170.    matFunc func( mName ) ; y.debugInfo( func ) ;
  171.    INDEX i, j, k, n = x.nCols(), nc = y.nCols() ;
  172.    DOUBLE r ;
  173.    error = NOERROR ;
  174.    if ( n > x.nRows() || n > y.nRows() ) {
  175.       error = GTDIM ;
  176.       return 0 ;
  177.    } // if
  178.    for ( k = 1 ; k <= nc ; k++ ) {
  179.       for ( i = n ; i ; i-- ) {
  180.          r = (DOUBLE) y(i,k) ;
  181.          // upper tri of x is upper tri of R
  182.          for ( j = i+1 ; j <= n ; j++ )
  183.             r -= x(i,j) * y(j,k) ;
  184.          if ( fabs( double( x(i,i) ) ) < tol ) {
  185.             error = ZEROD ;
  186.             return 0 ;
  187.          } // if
  188.          y(i,k) = (REAL) ( r / x(i,i) ) ;
  189.       } // for j
  190.    } // for k
  191.    return OK ;
  192. } // backSub
  193.  
  194. INDEX backSubT( matrix& y, const matrix& x, REAL tol,
  195.             matError& error )
  196. /*********************************************
  197.    Solve x'b = y assuming x is upper triang
  198.    Assume number of equations equals number
  199.    of columns in x. This allows x to have
  200.    more rows than columns e.g. after QR
  201.    decomposition. Proc requires the number
  202.    of rows be at least as great as the
  203.    number of equations.
  204.    Essentially the procedure is to use
  205.    forward substitution taking into account
  206.    x is not in transposed form.
  207. *********************************************/
  208. {
  209.    static char *mName = "backSubT" ;
  210.    matFunc func( mName ) ; y.debugInfo( func ) ;
  211.    INDEX i, j, k, n = x.nCols(), nc = y.nCols() ;
  212.    DOUBLE r ;
  213.    error = NOERROR ;
  214.    if ( n > x.nRows() || n > y.nRows() ) {
  215.       error = GTDIM ;
  216.       return 0 ;
  217.    } // if
  218.    for ( k = 1 ; k <= nc ; k++ ) {
  219.       for ( i = 1 ; i <= n ; i++ ) {
  220.          r = (DOUBLE) y(i,k) ;
  221.          // upper tri of x is upper tri of R
  222.          for ( j = 1 ; j < i ; j++ )
  223.             r -= x(j,i) * y(j,k) ;
  224.          if ( fabs( double( x(i,i) ) ) < tol ) {
  225.             error = ZEROD ;
  226.             return 0 ;
  227.          } // if
  228.          y(i,k) = (REAL) ( r / x(i,i) ) ;
  229.       } // for j
  230.    } // for k
  231.    return OK ;
  232. } // backSubT
  233.  
  234. INDEX triuInv( const matrix& R, matrix& Rinv, REAL tol,
  235.            matError& error )
  236. {
  237.    static char *mName = "triuInv" ;
  238.    matFunc func( mName ) ; R.debugInfo( func ) ;
  239.    INDEX i, j, k, n = R.nCols() ;
  240.    REAL sum ;
  241.    if ( n > R.nRows() ) {
  242.       error = GTDIM ;
  243.       return 0 ;
  244.    } // if
  245.    Rinv.reset( n, n ) ;
  246.    for ( j = 1 ; j <= n ; j++ ) {
  247.       if ( fabs( (double) R(j,j) ) < tol ) {
  248.          error = SINGM ;
  249.          return 0 ;
  250.       } // if
  251.       Rinv(j,j) = 1.0 / R(j,j) ;
  252.    } // for j
  253.    for ( j = n ; j ; j-- ) {
  254.       for ( i = j - 1 ; i ; i-- ) {
  255.          // find Rinv(i,j)
  256.          sum = 0.0 ;
  257.          for ( k = i + 1 ; k <= j ; k++ )
  258.             sum -= R(i,k) * Rinv(k,j) ;
  259.          Rinv(i,j) = Rinv(i,i) * sum ;
  260.       } // for i
  261.       for ( i = j + 1 ; i <= n ; i++ )
  262.          Rinv(i,j) = 0.0 ;
  263.    } // for j
  264.    return OK ;
  265. } // triuInv
  266.  
  267. INDEX qReflect( matrix& y, const matrix& q, matrix& w,
  268.                 matError& error )
  269. /******************************************************
  270.    Assume q and w are output from a Householder
  271.    QR decomposition of a matrix e.g. qrh. Now apply
  272.    reflectors in q to y. Insists on the number of rows
  273.    in q and y be equal. It also checks q and w have
  274.    the same number of rows. Note it handles the case
  275.    where the number of rows is less than or equal to
  276.    the number of columns.
  277. ******************************************************/
  278. {
  279.    static char *mName = "qReflect" ;
  280.    matFunc func( mName ) ; y.debugInfo( func ) ;
  281.    INDEX nr = q.nRows(), qnc = q.nCols(), ync = y.nCols() ;
  282.    refMatrix qCol(q), ySub(y) ;
  283.    if ( ( nr != y.nRows() ) || ( qnc != w.nRows() ) ) {
  284.       error = NEDIM ;
  285.       return 0 ;
  286.    } // if
  287.    INDEX k, n = ( nr <= qnc ) ? nr - 1 : qnc ;
  288.    for ( k = 1 ; k <= n ; k++ ) {
  289.       qCol.refSub( k, nr, k, k ) ;
  290.       ySub.refSub( k, nr, 1, ync ) ;
  291.       if ( !hReflect( ySub, qCol, w(k), error ) )
  292.          return 0 ;
  293.    } // for k
  294.    return OK ;
  295. } // qReflect
  296.  
  297. INDEX qTReflect( matrix& y, const matrix& q, matrix& w,
  298.          matError& error )
  299. /******************************************************
  300.    Assume q and w are output from a Householder
  301.    QR decomposition of a matrix e.g. qrh. Now apply
  302.    reflectors in q to y ***in reverse order***.
  303.    This is the inverse of qReflect.
  304.    Insists on q being square and the number of rows
  305.    in q and y be equal. It also checks q and w have
  306.    the same number of rows.
  307. ******************************************************/
  308. {
  309.    static char *mName = "qTReflect" ;
  310.    matFunc func( mName ) ; y.debugInfo( func ) ;
  311.    INDEX nr = q.nRows(), qnc = q.nCols(), ync = y.nCols() ;
  312.    refMatrix qCol(q), ySub(y) ;
  313.    if ( nr != qnc ) {
  314.       error = NTSQR ;
  315.       return 0 ;
  316.    } // if
  317.    if ( ( nr != y.nRows() ) || ( qnc != w.nRows() ) ) {
  318.       error = NEDIM ;
  319.       return 0 ;
  320.    } // if
  321.    INDEX k ;
  322.    for ( k = nr - 1 ; k ; k-- ) {
  323.       qCol.refSub( k, nr, k, k ) ;
  324.       ySub.refSub( k, nr, 1, ync ) ;
  325.       if ( !hReflect( ySub, qCol, w(k), error ) )
  326.          return 0 ;
  327.    } // for k
  328.    return OK ;
  329. } // qTReflect
  330.  
  331. INDEX hProject( matrix& y, matrix& x, matrix& R, matrix& w,
  332.                 matrix& beta, REAL tol, matError& error )
  333. /*************************************************************
  334.      Project y onto span(x). Use Householder QR of x.
  335.      Overwrite x with Q in partition form, R is returned
  336.      explicitly. w has the scales for the reflectors i.e.
  337.      2/u'u. y overwritten by projected y, beta with OLS
  338.      coeffs.
  339. *************************************************************/
  340. {
  341.    static char *mName = "hProject" ;
  342.    matFunc func( mName ) ; y.debugInfo( func ) ;
  343.    INDEX nr = x.nRows(), xnc = x.nCols(), ync = y.nCols() ;
  344.    matrix diag ;
  345.    error = NOERROR ;
  346.    if ( nr != y.nRows() ) {
  347.       error = NEDIM ;
  348.       return 0 ;
  349.    } // if
  350.    // form QR decomposition of x
  351.    // Q returned in lower x and w, R is upper x and diag
  352.    if ( !qrh( x, diag, w, tol, error ) )
  353.       return 0 ;
  354.    // Apply householder transforms to y
  355.    if ( !qReflect( y, x, w, error ) )
  356.       return 0 ;
  357.    // form R explicitly
  358.    INDEX n = ( nr <= xnc ) ? nr - 1 : xnc ;
  359.    R.triuOf( x.sub( 1, n, 1, n ) ) ;
  360.    R.setDiag( diag ) ;
  361.    beta = y.sub( 1, n, 1, ync ) ;
  362.    // Apply backSub to solve Rb = beta.
  363.    // leaving OLS coeffs in upper beta
  364.    if ( !backSub( beta, R, tol, error ) )
  365.       return 0 ;
  366.    return OK ;
  367. } // hProject
  368.