home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / mfuntort < prev    next >
Encoding:
Text File  |  1994-01-13  |  4.4 KB  |  182 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. /* mfuntort.c,  10/11/93 */
  28.  
  29. static char rcsid[] = "$Id: mfuntort.c,v 1.2 1994/01/14 01:08:06 des Exp $";
  30.  
  31. #include        <stdio.h>
  32. #include        <math.h>
  33. #include        "matrix.h"
  34. #include        "matrix2.h"
  35.  
  36.  
  37. #define errmesg(mesg)   printf("Error: %s error: line %d\n",mesg,__LINE__)
  38. #define notice(mesg)    printf("# Testing %s...\n",mesg);
  39.  
  40. #define DIM  10
  41.  
  42. void main()
  43. {
  44.  
  45.    MAT *A, *B, *C, *OUTA, *OUTB, *TMP;
  46.    MAT *exp_A_expected, *exp_A;
  47.    VEC *x, *b;
  48.    double c, eps = 1e-10;
  49.    int i, j, q_out, j_out;
  50.  
  51.    mem_info_on(TRUE);
  52.  
  53.    A = m_get(DIM,DIM);
  54.    B = m_get(DIM,DIM);
  55.    C = m_get(DIM,DIM);
  56.    OUTA = m_get(DIM,DIM);
  57.    OUTB = m_get(DIM,DIM);
  58.    TMP = m_get(DIM,DIM);
  59.    x = v_get(DIM);
  60.    b = v_get(6);
  61.  
  62.    notice("exponent of a matrix");
  63.  
  64.    m_ident(A);
  65.    mem_stat_mark(1);
  66.    _m_exp(A,eps,OUTA,&q_out,&j_out);
  67.    printf("# q_out = %d, j_out = %d\n",q_out,j_out);
  68.  
  69.    m_exp(A,eps,OUTA);
  70.    sm_mlt(exp(1.0),A,A);
  71.    m_sub(OUTA,A,TMP);
  72.    printf("# ||exp(I) - e*I|| = %g\n",m_norm_inf(TMP));
  73.  
  74.    m_rand(A);
  75.    m_transp(A,TMP);
  76.    m_add(A,TMP,A);
  77.    B = m_copy(A,B);
  78.  
  79.    m_exp(A,eps,OUTA);
  80.  
  81.    symmeig(B,OUTB,x);
  82.    m_zero(TMP);
  83.    for (i=0; i < x->dim; i++)
  84.      TMP->me[i][i] = exp(x->ve[i]);
  85.    m_mlt(OUTB,TMP,C);
  86.    mmtr_mlt(C,OUTB,TMP);
  87.    m_sub(TMP,OUTA,TMP);
  88.    printf("# ||exp(A) - Q*exp(lambda)*Q^T|| = %g\n",m_norm_inf(TMP));
  89.  
  90.    notice("polynomial of a matrix");
  91.    m_rand(A);
  92.    m_transp(A,TMP);
  93.    m_add(A,TMP,A);
  94.    B = m_copy(A,B);
  95.    v_rand(b);
  96.  
  97.    m_poly(A,b,OUTA);
  98.  
  99.    symmeig(B,OUTB,x);
  100.    m_zero(TMP);
  101.    for (i=0; i < x->dim; i++) {
  102.       c = b->ve[b->dim-1];
  103.       for (j=b->dim-2; j >= 0; j--) 
  104.     c = c*x->ve[i] + b->ve[j];
  105.       TMP->me[i][i] = c;
  106.    }
  107.    m_mlt(OUTB,TMP,C);
  108.    mmtr_mlt(C,OUTB,TMP);
  109.    m_sub(TMP,OUTA,TMP);
  110.    printf("# ||poly(A) - Q*poly(lambda)*Q^T|| = %g\n",m_norm_inf(TMP));
  111.    mem_stat_free(1);
  112.  
  113.  
  114.    /* Brook Milligan's test */
  115.  
  116.    M_FREE(A);
  117.    M_FREE(B);
  118.    M_FREE(C);
  119.  
  120.    notice("exponent of a nonsymmetric matrix");
  121.    A = m_get (2, 2);
  122.    A -> me [0][0] = 1.0;
  123.    A -> me [0][1] = 1.0;
  124.    A -> me [1][0] = 4.0;
  125.    A -> me [1][1] = 1.0;
  126.    
  127.    exp_A_expected = m_get(2, 2);
  128.    exp_A_expected -> me [0][0] = exp (3.0) / 2.0 + exp (-1.0) / 2.0;
  129.    exp_A_expected -> me [0][1] = exp (3.0) / 4.0 - exp (-1.0) / 4.0;
  130.    exp_A_expected -> me [1][0] = exp (3.0)       - exp (-1.0);
  131.    exp_A_expected -> me [1][1] = exp (3.0) / 2.0 + exp (-1.0) / 2.0;
  132.    
  133.    printf ("A:\n");
  134.    for (i = 0; i < 2; i++)
  135.    {
  136.       for (j = 0; j < 2; j++)
  137.         printf ("   %15.8e", A -> me [i][j]);
  138.       printf ("\n");
  139.    }
  140.    
  141.    printf ("\nexp(A) (expected):\n");
  142.    for (i = 0; i < 2; i++)
  143.    {
  144.       for (j = 0; j < 2; j++)
  145.         printf ("   %15.8e", exp_A_expected -> me [i][j]);
  146.       printf ("\n");
  147.    }
  148.    
  149.    mem_stat_mark(3);
  150.    exp_A = m_exp (A, 1e-16,NULL);
  151.    mem_stat_free(3);
  152.  
  153.    printf ("\nexp(A):\n");
  154.    for (i = 0; i < 2; i++)
  155.    {
  156.       for (j = 0; j < 2; j++)
  157.         printf ("   %15.8e", exp_A -> me [i][j]);
  158.       printf ("\n");
  159.    }
  160.    printf ("\nexp(A) - exp(A) (expected):\n");
  161.    for (i = 0; i < 2; i++)
  162.    {
  163.       for (j = 0; j < 2; j++)
  164.         printf ("   %15.8e", exp_A -> me [i][j] - exp_A_expected -> me [i][j]);
  165.       printf ("\n");
  166.    }
  167.  
  168.    M_FREE(A);
  169.    M_FREE(B);
  170.    M_FREE(C);
  171.    M_FREE(exp_A);
  172.    M_FREE(exp_A_expected);
  173.    M_FREE(OUTA);
  174.    M_FREE(OUTB);
  175.    M_FREE(TMP);
  176.    V_FREE(b);
  177.    V_FREE(x);
  178.  
  179.    mem_info();
  180. }
  181.  
  182.