home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / zlufctr < prev    next >
Encoding:
Text File  |  1994-01-12  |  6.7 KB  |  279 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.     Matrix factorisation routines to work with the other matrix files.
  28.     Complex version
  29. */
  30.  
  31. static    char    rcsid[] = "$Id: zlufctr.c,v 1.1 1994/01/13 04:26:20 des Exp $";
  32.  
  33. #include    <stdio.h>
  34. #include    <math.h>
  35. #include    "zmatrix.h"
  36. #include        "zmatrix2.h"
  37.  
  38. #define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)
  39.  
  40.  
  41. /* Most matrix factorisation routines are in-situ unless otherwise specified */
  42.  
  43. /* zLUfactor -- Gaussian elimination with scaled partial pivoting
  44.         -- Note: returns LU matrix which is A */
  45. ZMAT    *zLUfactor(A,pivot)
  46. ZMAT    *A;
  47. PERM    *pivot;
  48. {
  49.     u_int    i, j, k, k_max, m, n;
  50.     int    i_max;
  51.     Real    dtemp, max1;
  52.     complex    **A_v, *A_piv, *A_row, temp;
  53.     static    VEC    *scale = VNULL;
  54.  
  55.     if ( A==ZMNULL || pivot==PNULL )
  56.         error(E_NULL,"zLUfactor");
  57.     if ( pivot->size != A->m )
  58.         error(E_SIZES,"zLUfactor");
  59.     m = A->m;    n = A->n;
  60.     scale = v_resize(scale,A->m);
  61.     MEM_STAT_REG(scale,TYPE_VEC);
  62.     A_v = A->me;
  63.  
  64.     /* initialise pivot with identity permutation */
  65.     for ( i=0; i<m; i++ )
  66.         pivot->pe[i] = i;
  67.  
  68.     /* set scale parameters */
  69.     for ( i=0; i<m; i++ )
  70.     {
  71.         max1 = 0.0;
  72.         for ( j=0; j<n; j++ )
  73.         {
  74.             dtemp = zabs(A_v[i][j]);
  75.             max1 = max(max1,dtemp);
  76.         }
  77.         scale->ve[i] = max1;
  78.     }
  79.  
  80.     /* main loop */
  81.     k_max = min(m,n)-1;
  82.     for ( k=0; k<k_max; k++ )
  83.     {
  84.         /* find best pivot row */
  85.         max1 = 0.0;    i_max = -1;
  86.         for ( i=k; i<m; i++ )
  87.         if ( scale->ve[i] > 0.0 )
  88.         {
  89.             dtemp = zabs(A_v[i][k])/scale->ve[i];
  90.             if ( dtemp > max1 )
  91.             { max1 = dtemp;    i_max = i;    }
  92.         }
  93.         
  94.         /* if no pivot then ignore column k... */
  95.         if ( i_max == -1 )
  96.         continue;
  97.  
  98.         /* do we pivot ? */
  99.         if ( i_max != k )    /* yes we do... */
  100.         {
  101.         px_transp(pivot,i_max,k);
  102.         for ( j=0; j<n; j++ )
  103.         {
  104.             temp = A_v[i_max][j];
  105.             A_v[i_max][j] = A_v[k][j];
  106.             A_v[k][j] = temp;
  107.         }
  108.         }
  109.         
  110.         /* row operations */
  111.         for ( i=k+1; i<m; i++ )    /* for each row do... */
  112.         {    /* Note: divide by zero should never happen */
  113.         temp = A_v[i][k] = zdiv(A_v[i][k],A_v[k][k]);
  114.         A_piv = &(A_v[k][k+1]);
  115.         A_row = &(A_v[i][k+1]);
  116.         temp.re = - temp.re;
  117.         temp.im = - temp.im;
  118.         if ( k+1 < n )
  119.             __zmltadd__(A_row,A_piv,temp,(int)(n-(k+1)),Z_NOCONJ);
  120.         /*********************************************
  121.           for ( j=k+1; j<n; j++ )
  122.           A_v[i][j] -= temp*A_v[k][j];
  123.           (*A_row++) -= temp*(*A_piv++);
  124.         *********************************************/
  125.         }
  126.     }
  127.  
  128.     return A;
  129. }
  130.  
  131.  
  132. /* zLUsolve -- given an LU factorisation in A, solve Ax=b */
  133. ZVEC    *zLUsolve(A,pivot,b,x)
  134. ZMAT    *A;
  135. PERM    *pivot;
  136. ZVEC    *b,*x;
  137. {
  138.     if ( A==ZMNULL || b==ZVNULL || pivot==PNULL )
  139.         error(E_NULL,"zLUsolve");
  140.     if ( A->m != A->n || A->n != b->dim )
  141.         error(E_SIZES,"zLUsolve");
  142.  
  143.     x = px_zvec(pivot,b,x);    /* x := P.b */
  144.     zLsolve(A,x,x,1.0);    /* implicit diagonal = 1 */
  145.     zUsolve(A,x,x,0.0);    /* explicit diagonal */
  146.  
  147.     return (x);
  148. }
  149.  
  150. /* zLUAsolve -- given an LU factorisation in A, solve A^*.x=b */
  151. ZVEC    *zLUAsolve(LU,pivot,b,x)
  152. ZMAT    *LU;
  153. PERM    *pivot;
  154. ZVEC    *b,*x;
  155. {
  156.     if ( ! LU || ! b || ! pivot )
  157.         error(E_NULL,"zLUAsolve");
  158.     if ( LU->m != LU->n || LU->n != b->dim )
  159.         error(E_SIZES,"zLUAsolve");
  160.  
  161.     x = zv_copy(b,x);
  162.     zUAsolve(LU,x,x,0.0);    /* explicit diagonal */
  163.     zLAsolve(LU,x,x,1.0);    /* implicit diagonal = 1 */
  164.     pxinv_zvec(pivot,x,x);    /* x := P^*.x */
  165.  
  166.     return (x);
  167. }
  168.  
  169. /* zm_inverse -- returns inverse of A, provided A is not too rank deficient
  170.     -- uses LU factorisation */
  171. ZMAT    *zm_inverse(A,out)
  172. ZMAT    *A, *out;
  173. {
  174.     int    i;
  175.     ZVEC    *tmp, *tmp2;
  176.     ZMAT    *A_cp;
  177.     PERM    *pivot;
  178.  
  179.     if ( ! A )
  180.         error(E_NULL,"zm_inverse");
  181.     if ( A->m != A->n )
  182.         error(E_SQUARE,"zm_inverse");
  183.     if ( ! out || out->m < A->m || out->n < A->n )
  184.         out = zm_resize(out,A->m,A->n);
  185.  
  186.     A_cp = zm_copy(A,ZMNULL);
  187.     tmp = zv_get(A->m);
  188.     tmp2 = zv_get(A->m);
  189.     pivot = px_get(A->m);
  190.     tracecatch(zLUfactor(A_cp,pivot),"zm_inverse");
  191.     for ( i = 0; i < A->n; i++ )
  192.     {
  193.         zv_zero(tmp);
  194.         tmp->ve[i].re = 1.0;
  195.         tmp->ve[i].im = 0.0;
  196.         tracecatch(zLUsolve(A_cp,pivot,tmp,tmp2),"m_inverse");
  197.         zset_col(out,i,tmp2);
  198.     }
  199.  
  200.     ZM_FREE(A_cp);
  201.     ZV_FREE(tmp);    ZV_FREE(tmp2);
  202.     PX_FREE(pivot);
  203.  
  204.     return out;
  205. }
  206.  
  207. /* zLUcondest -- returns an estimate of the condition number of LU given the
  208.     LU factorisation in compact form */
  209. double    zLUcondest(LU,pivot)
  210. ZMAT    *LU;
  211. PERM    *pivot;
  212. {
  213.     static    ZVEC    *y = ZVNULL, *z = ZVNULL;
  214.     Real    cond_est, L_norm, U_norm, norm, sn_inv;
  215.     complex    sum;
  216.     int        i, j, n;
  217.  
  218.     if ( ! LU || ! pivot )
  219.     error(E_NULL,"zLUcondest");
  220.     if ( LU->m != LU->n )
  221.     error(E_SQUARE,"zLUcondest");
  222.     if ( LU->n != pivot->size )
  223.     error(E_SIZES,"zLUcondest");
  224.  
  225.     n = LU->n;
  226.     y = zv_resize(y,n);
  227.     z = zv_resize(z,n);
  228.     MEM_STAT_REG(y,TYPE_ZVEC);
  229.     MEM_STAT_REG(z,TYPE_ZVEC);
  230.  
  231.     cond_est = 0.0;        /* should never be returned */
  232.  
  233.     for ( i = 0; i < n; i++ )
  234.     {
  235.     sum.re = 1.0;
  236.     sum.im = 0.0;
  237.     for ( j = 0; j < i; j++ )
  238.         /* sum -= LU->me[j][i]*y->ve[j]; */
  239.         sum = zsub(sum,zmlt(LU->me[j][i],y->ve[j]));
  240.     /* sum -= (sum < 0.0) ? 1.0 : -1.0; */
  241.     sn_inv = 1.0 / zabs(sum);
  242.     sum.re += sum.re * sn_inv;
  243.     sum.im += sum.im * sn_inv;
  244.     if ( is_zero(LU->me[i][i]) )
  245.         return HUGE;
  246.     /* y->ve[i] = sum / LU->me[i][i]; */
  247.     y->ve[i] = zdiv(sum,LU->me[i][i]);
  248.     }
  249.  
  250.     zLAsolve(LU,y,y,1.0);
  251.     zLUsolve(LU,pivot,y,z);
  252.  
  253.     /* now estimate norm of A (even though it is not directly available) */
  254.     /* actually computes ||L||_inf.||U||_inf */
  255.     U_norm = 0.0;
  256.     for ( i = 0; i < n; i++ )
  257.     {
  258.     norm = 0.0;
  259.     for ( j = i; j < n; j++ )
  260.         norm += zabs(LU->me[i][j]);
  261.     if ( norm > U_norm )
  262.         U_norm = norm;
  263.     }
  264.     L_norm = 0.0;
  265.     for ( i = 0; i < n; i++ )
  266.     {
  267.     norm = 1.0;
  268.     for ( j = 0; j < i; j++ )
  269.         norm += zabs(LU->me[i][j]);
  270.     if ( norm > L_norm )
  271.         L_norm = norm;
  272.     }
  273.  
  274.     tracecatch(cond_est = U_norm*L_norm*zv_norm_inf(z)/zv_norm_inf(y),
  275.            "LUcondest");
  276.  
  277.     return cond_est;
  278. }
  279.