home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / zsolve.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  8KB  |  301 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.     Complex case
  30. */
  31.  
  32. static    char    rcsid[] = "$Id: zsolve.c,v 1.1 1994/01/13 04:20:33 des Exp $";
  33.  
  34. #include    <stdio.h>
  35. #include    <math.h>
  36. #include        "zmatrix2.h"
  37.  
  38.  
  39. #define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0 )
  40.  
  41. /* Most matrix factorisation routines are in-situ unless otherwise specified */
  42.  
  43. /* zUsolve -- back substitution with optional over-riding diagonal
  44.         -- can be in-situ but doesn't need to be */
  45. ZVEC    *zUsolve(matrix,b,out,diag)
  46. ZMAT    *matrix;
  47. ZVEC    *b, *out;
  48. double    diag;
  49. {
  50.     u_int    dim /* , j */;
  51.     int    i, i_lim;
  52.     complex    **mat_ent, *mat_row, *b_ent, *out_ent, *out_col, sum;
  53.     
  54.     if ( matrix==ZMNULL || b==ZVNULL )
  55.     error(E_NULL,"zUsolve");
  56.     dim = min(matrix->m,matrix->n);
  57.     if ( b->dim < dim )
  58.     error(E_SIZES,"zUsolve");
  59.     if ( out==ZVNULL || out->dim < dim )
  60.     out = zv_resize(out,matrix->n);
  61.     mat_ent = matrix->me;    b_ent = b->ve;    out_ent = out->ve;
  62.     
  63.     for ( i=dim-1; i>=0; i-- )
  64.     if ( ! is_zero(b_ent[i]) )
  65.         break;
  66.     else
  67.         out_ent[i].re = out_ent[i].im = 0.0;
  68.     i_lim = i;
  69.     
  70.     for ( i = i_lim; i>=0; i-- )
  71.     {
  72.     sum = b_ent[i];
  73.     mat_row = &(mat_ent[i][i+1]);
  74.     out_col = &(out_ent[i+1]);
  75.     sum = zsub(sum,__zip__(mat_row,out_col,i_lim-i,Z_NOCONJ));
  76.     /******************************************************
  77.       for ( j=i+1; j<=i_lim; j++ )
  78.       sum -= mat_ent[i][j]*out_ent[j];
  79.       sum -= (*mat_row++)*(*out_col++);
  80.     ******************************************************/
  81.     if ( diag == 0.0 )
  82.     {
  83.         if ( is_zero(mat_ent[i][i]) )
  84.         error(E_SING,"zUsolve");
  85.         else
  86.         /* out_ent[i] = sum/mat_ent[i][i]; */
  87.         out_ent[i] = zdiv(sum,mat_ent[i][i]);
  88.     }
  89.     else
  90.     {
  91.         /* out_ent[i] = sum/diag; */
  92.         out_ent[i].re = sum.re / diag;
  93.         out_ent[i].im = sum.im / diag;
  94.     }
  95.     }
  96.     
  97.     return (out);
  98. }
  99.  
  100. /* zLsolve -- forward elimination with (optional) default diagonal value */
  101. ZVEC    *zLsolve(matrix,b,out,diag)
  102. ZMAT    *matrix;
  103. ZVEC    *b,*out;
  104. double    diag;
  105. {
  106.     u_int    dim, i, i_lim /* , j */;
  107.     complex    **mat_ent, *mat_row, *b_ent, *out_ent, *out_col, sum;
  108.     
  109.     if ( matrix==ZMNULL || b==ZVNULL )
  110.     error(E_NULL,"zLsolve");
  111.     dim = min(matrix->m,matrix->n);
  112.     if ( b->dim < dim )
  113.     error(E_SIZES,"zLsolve");
  114.     if ( out==ZVNULL || out->dim < dim )
  115.     out = zv_resize(out,matrix->n);
  116.     mat_ent = matrix->me;    b_ent = b->ve;    out_ent = out->ve;
  117.     
  118.     for ( i=0; i<dim; i++ )
  119.     if ( ! is_zero(b_ent[i]) )
  120.         break;
  121.     else
  122.         out_ent[i].re = out_ent[i].im = 0.0;
  123.     i_lim = i;
  124.     
  125.     for ( i = i_lim; i<dim; i++ )
  126.     {
  127.     sum = b_ent[i];
  128.     mat_row = &(mat_ent[i][i_lim]);
  129.     out_col = &(out_ent[i_lim]);
  130.     sum = zsub(sum,__zip__(mat_row,out_col,(int)(i-i_lim),Z_NOCONJ));
  131.     /*****************************************************
  132.       for ( j=i_lim; j<i; j++ )
  133.       sum -= mat_ent[i][j]*out_ent[j];
  134.       sum -= (*mat_row++)*(*out_col++);
  135.     ******************************************************/
  136.     if ( diag == 0.0 )
  137.     {
  138.         if ( is_zero(mat_ent[i][i]) )
  139.         error(E_SING,"zLsolve");
  140.         else
  141.         out_ent[i] = zdiv(sum,mat_ent[i][i]);
  142.     }
  143.     else
  144.     {
  145.         out_ent[i].re = sum.re / diag;
  146.         out_ent[i].im = sum.im / diag;
  147.     }
  148.     }
  149.     
  150.     return (out);
  151. }
  152.  
  153.  
  154. /* zUAsolve -- forward elimination with (optional) default diagonal value
  155.         using UPPER triangular part of matrix */
  156. ZVEC    *zUAsolve(U,b,out,diag)
  157. ZMAT    *U;
  158. ZVEC    *b,*out;
  159. double    diag;
  160. {
  161.     u_int    dim, i, i_lim /* , j */;
  162.     complex    **U_me, *b_ve, *out_ve, tmp;
  163.     Real    invdiag;
  164.     
  165.     if ( ! U || ! b )
  166.     error(E_NULL,"zUAsolve");
  167.     dim = min(U->m,U->n);
  168.     if ( b->dim < dim )
  169.     error(E_SIZES,"zUAsolve");
  170.     out = zv_resize(out,U->n);
  171.     U_me = U->me;    b_ve = b->ve;    out_ve = out->ve;
  172.     
  173.     for ( i=0; i<dim; i++ )
  174.     if ( ! is_zero(b_ve[i]) )
  175.         break;
  176.     else
  177.         out_ve[i].re = out_ve[i].im = 0.0;
  178.     i_lim = i;
  179.     if ( b != out )
  180.     {
  181.     __zzero__(out_ve,out->dim);
  182.     /* MEM_COPY(&(b_ve[i_lim]),&(out_ve[i_lim]),
  183.        (dim-i_lim)*sizeof(complex)); */
  184.     MEMCOPY(&(b_ve[i_lim]),&(out_ve[i_lim]),dim-i_lim,complex);
  185.     }
  186.  
  187.     if ( diag == 0.0 )
  188.     {
  189.     for (    ; i<dim; i++ )
  190.     {
  191.         tmp = zconj(U_me[i][i]);
  192.         if ( is_zero(tmp) )
  193.         error(E_SING,"zUAsolve");
  194.         /* out_ve[i] /= tmp; */
  195.         out_ve[i] = zdiv(out_ve[i],tmp);
  196.         tmp.re = - out_ve[i].re;
  197.         tmp.im = - out_ve[i].im;
  198.         __zmltadd__(&(out_ve[i+1]),&(U_me[i][i+1]),tmp,dim-i-1,Z_CONJ);
  199.     }
  200.     }
  201.     else
  202.     {
  203.     invdiag = 1.0/diag;
  204.     for (    ; i<dim; i++ )
  205.     {
  206.         out_ve[i].re *= invdiag;
  207.         out_ve[i].im *= invdiag;
  208.         tmp.re = - out_ve[i].re;
  209.         tmp.im = - out_ve[i].im;
  210.         __zmltadd__(&(out_ve[i+1]),&(U_me[i][i+1]),tmp,dim-i-1,Z_CONJ);
  211.     }
  212.     }
  213.     return (out);
  214. }
  215.  
  216. /* zDsolve -- solves Dx=b where D is the diagonal of A -- may be in-situ */
  217. ZVEC    *zDsolve(A,b,x)
  218. ZMAT    *A;
  219. ZVEC    *b,*x;
  220. {
  221.     u_int    dim, i;
  222.     
  223.     if ( ! A || ! b )
  224.     error(E_NULL,"zDsolve");
  225.     dim = min(A->m,A->n);
  226.     if ( b->dim < dim )
  227.     error(E_SIZES,"zDsolve");
  228.     x = zv_resize(x,A->n);
  229.     
  230.     dim = b->dim;
  231.     for ( i=0; i<dim; i++ )
  232.     if ( is_zero(A->me[i][i]) )
  233.         error(E_SING,"zDsolve");
  234.     else
  235.         x->ve[i] = zdiv(b->ve[i],A->me[i][i]);
  236.     
  237.     return (x);
  238. }
  239.  
  240. /* zLAsolve -- back substitution with optional over-riding diagonal
  241.         using the LOWER triangular part of matrix
  242.         -- can be in-situ but doesn't need to be */
  243. ZVEC    *zLAsolve(L,b,out,diag)
  244. ZMAT    *L;
  245. ZVEC    *b, *out;
  246. double    diag;
  247. {
  248.     u_int    dim;
  249.     int        i, i_lim;
  250.     complex    **L_me, *b_ve, *out_ve, tmp;
  251.     Real    invdiag;
  252.     
  253.     if ( ! L || ! b )
  254.     error(E_NULL,"zLAsolve");
  255.     dim = min(L->m,L->n);
  256.     if ( b->dim < dim )
  257.     error(E_SIZES,"zLAsolve");
  258.     out = zv_resize(out,L->n);
  259.     L_me = L->me;    b_ve = b->ve;    out_ve = out->ve;
  260.     
  261.     for ( i=dim-1; i>=0; i-- )
  262.     if ( ! is_zero(b_ve[i]) )
  263.         break;
  264.     i_lim = i;
  265.  
  266.     if ( b != out )
  267.     {
  268.     __zzero__(out_ve,out->dim);
  269.     /* MEM_COPY(b_ve,out_ve,(i_lim+1)*sizeof(complex)); */
  270.     MEMCOPY(b_ve,out_ve,i_lim+1,complex);
  271.     }
  272.  
  273.     if ( diag == 0.0 )
  274.     {
  275.     for (        ; i>=0; i-- )
  276.     {
  277.         tmp = zconj(L_me[i][i]);
  278.         if ( is_zero(tmp) )
  279.         error(E_SING,"zLAsolve");
  280.         out_ve[i] = zdiv(out_ve[i],tmp);
  281.         tmp.re = - out_ve[i].re;
  282.         tmp.im = - out_ve[i].im;
  283.         __zmltadd__(out_ve,L_me[i],tmp,i,Z_CONJ);
  284.     }
  285.     }
  286.     else
  287.     {
  288.     invdiag = 1.0/diag;
  289.     for (        ; i>=0; i-- )
  290.     {
  291.         out_ve[i].re *= invdiag;
  292.         out_ve[i].im *= invdiag;
  293.         tmp.re = - out_ve[i].re;
  294.         tmp.im = - out_ve[i].im;
  295.         __zmltadd__(out_ve,L_me[i],tmp,i,Z_CONJ);
  296.     }
  297.     }
  298.     
  299.     return (out);
  300. }
  301.