home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / update < prev    next >
Encoding:
Text File  |  1994-01-13  |  3.4 KB  |  132 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*
  28.     Matrix factorisation routines to work with the other matrix files.
  29. */
  30.  
  31. /* update.c 1.3 11/25/87 */
  32. static    char    rcsid[] = "$Id: update.c,v 1.2 1994/01/13 05:26:06 des Exp $";
  33.  
  34. #include    <stdio.h>
  35. #include    <math.h>
  36. #include    "matrix.h"
  37. #include        "matrix2.h"
  38.  
  39.  
  40.  
  41.  
  42. /* Most matrix factorisation routines are in-situ unless otherwise specified */
  43.  
  44. /* LDLupdate -- updates a CHolesky factorisation, replacing LDL' by
  45.     MD~M' = LDL' + alpha.w.w' Note: w is overwritten
  46.     Ref: Gill et al Math Comp 28, p516 Algorithm C1 */
  47. MAT    *LDLupdate(CHmat,w,alpha)
  48. MAT    *CHmat;
  49. VEC    *w;
  50. double    alpha;
  51. {
  52.     u_int    i,j;
  53.     Real    diag,new_diag,beta,p;
  54.  
  55.     if ( CHmat==(MAT *)NULL || w==(VEC *)NULL )
  56.         error(E_NULL,"LDLupdate");
  57.     if ( CHmat->m != CHmat->n || w->dim != CHmat->m )
  58.         error(E_SIZES,"LDLupdate");
  59.  
  60.     for ( j=0; j < w->dim; j++ )
  61.     {
  62.         p = w->ve[j];
  63.         diag = CHmat->me[j][j];
  64.         new_diag = CHmat->me[j][j] = diag + alpha*p*p;
  65.         if ( new_diag <= 0.0 )
  66.             error(E_POSDEF,"LDLupdate");
  67.         beta = p*alpha/new_diag;
  68.         alpha *= diag/new_diag;
  69.  
  70.         for ( i=j+1; i < w->dim; i++ )
  71.         {
  72.             w->ve[i] -= p*CHmat->me[i][j];
  73.             CHmat->me[i][j] += beta*w->ve[i];
  74.             CHmat->me[j][i] = CHmat->me[i][j];
  75.         }
  76.     }
  77.  
  78.     return (CHmat);
  79. }
  80.  
  81.  
  82. /* QRupdate -- updates QR factorisation in expanded form (seperate matrices)
  83.     Finds Q+, R+ s.t. Q+.R+ = Q.(R+u.v') and Q+ orthogonal, R+ upper triang
  84.     Ref: Golub & van Loan Matrix Computations pp437-443
  85.     -- does not update Q if it is NULL */
  86. MAT    *QRupdate(Q,R,u,v)
  87. MAT    *Q,*R;
  88. VEC    *u,*v;
  89. {
  90.     int    i,j,k;
  91.     Real    c,s,temp;
  92.  
  93.     if ( ! R || ! u || ! v )
  94.         error(E_NULL,"QRupdate");
  95.     if ( ( Q && ( Q->m != Q->n || R->m != Q->n ) ) ||
  96.                     u->dim != R->m || v->dim != R->n )
  97.         error(E_SIZES,"QRupdate");
  98.  
  99.     /* find largest k s.t. u[k] != 0 */
  100.     for ( k=R->m-1; k>=0; k-- )
  101.         if ( u->ve[k] != 0.0 )
  102.             break;
  103.  
  104.     /* transform R+u.v' to Hessenberg form */
  105.     for ( i=k-1; i>=0; i-- )
  106.     {
  107.         /* get Givens rotation */
  108.         givens(u->ve[i],u->ve[i+1],&c,&s);
  109.         rot_rows(R,i,i+1,c,s,R);
  110.         if ( Q )
  111.             rot_cols(Q,i,i+1,c,s,Q);
  112.         rot_vec(u,i,i+1,c,s,u);
  113.     }
  114.  
  115.     /* add into R */
  116.     temp = u->ve[0];
  117.     for ( j=0; j<R->n; j++ )
  118.         R->me[0][j] += temp*v->ve[j];
  119.  
  120.     /* transform Hessenberg to upper triangular */
  121.     for ( i=0; i<k; i++ )
  122.     {
  123.         givens(R->me[i][i],R->me[i+1][i],&c,&s);
  124.         rot_rows(R,i,i+1,c,s,R);
  125.         if ( Q )
  126.             rot_cols(Q,i,i+1,c,s,Q);
  127.     }
  128.  
  129.     return R;
  130. }
  131.  
  132.