home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / zmatlab < prev    next >
Encoding:
Text File  |  1994-01-12  |  5.7 KB  |  222 lines
  1.    mtrm_mlt(Q,Q,D);
  2.     for ( i = 0; i < D->m; i++ )
  3.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  4.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*3 )
  5.     {
  6.     errmesg("symmeig()");
  7.     printf("# symmeig() orthogonality error = %g [cf MACHEPS = %g]\n",
  8.            m_norm1(D), MACHEPS);
  9.     }
  10.  
  11.     MEMCHK();
  12.  
  13.     /* now test (real) Schur decomposition */
  14.     /* m_copy(A,B); */
  15.     M_FREE(A);
  16.     A = m_get(11,11);
  17.     m_rand(A);
  18.     B = m_copy(A,B);
  19.     MEMCHK();
  20.  
  21.     B = schur(B,Q);
  22.     MEMCHK();
  23.  
  24.     m_mlt(Q,B,C);
  25.     mmtr_mlt(C,Q,D);
  26.     MEMCHK();
  27.     m_sub(A,D,D);
  28.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*m_norm1(B)*5 )
  29.     {
  30.     errmesg("schur()");
  31.     printf("# Schur reconstruction error = %g [cf MACHEPS = %g]\n",
  32.            m_norm1(D), MACHEPS);
  33.     }
  34.  
  35.     /* orthogonality check */
  36.     mmtr_mlt(Q,Q,D);
  37.     for ( i = 0; i < D->m; i++ )
  38.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  39.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*10 )
  40.     {
  41.     errmesg("schur()");
  42.     printf("# Schur orthogonality error = %g [cf MACHEPS = %g]\n",
  43.            m_norm1(D), MACHEPS);
  44.     }
  45.  
  46.     MEMCHK();
  47.  
  48.     /* check for termination */
  49.     Atmp = m_get(2,2);
  50.     Qtmp = m_get(2,2);
  51.     /* this is a 2 x 2 matrix with real eigenvalues */
  52.     Atmp->me[0][0] = 1;
  53.     Atmp->me[1][1] = 1;
  54.     Atmp->me[0][1] = 4;
  55.     Atmp->me[1][0] = 1;
  56.     schur(Atmp,Qtmp);
  57.  
  58.     MEMCHK();
  59.  
  60.     /* now test SVD */
  61.     A = m_resize(A,11,7);
  62.     m_rand(A);
  63.     U = m_get(A->n,A->n);
  64.     Q = m_resize(Q,A->m,A->m);
  65.     u = v_resize(u,max(A->m,A->n));
  66.     svd(A,Q,U,u);
  67.     /* check reconstruction of A */
  68.     D = m_resize(D,A->m,A->n);
  69.     C = m_resize(C,A->m,A->n);
  70.     m_zero(D);
  71.     for ( i = 0; i < min(A->m,A->n); i++ )
  72.     m_set_val(D,i,i,v_entry(u,i));
  73.     mtrm_mlt(Q,D,C);
  74.     m_mlt(C,U,D);
  75.     m_sub(A,D,D);
  76.     if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(Q)*m_norm1(A) )
  77.     {
  78.     errmesg("svd()");
  79.     printf("# SVD reconstruction error = %g [cf MACHEPS = %g]\n",
  80.            m_norm1(D), MACHEPS);
  81.     }
  82.     /* check orthogonality of Q and U */
  83.     D = m_resize(D,Q->n,Q->n);
  84.     mtrm_mlt(Q,Q,D);
  85.     for ( i = 0; i < D->m; i++ )
  86.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  87.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*5 )
  88.     {
  89.     errmesg("svd()");
  90.     printf("# SVD orthognality error (Q) = %g [cf MACHEPS = %g\n",
  91.            m_norm1(D), MACHEPS);
  92.     }
  93.     D = m_resize(D,U->n,U->n);
  94.     mtrm_mlt(U,U,D);
  95.     for ( i = 0; i < D->m; i++ )
  96.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  97.     if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(U)*5 )
  98.     {
  99.     errmesg("svd()");
  100.     printf("# SVD orthognality error (U) = %g [cf MACHEPS = %g\n",
  101.            m_norm1(D), MACHEPS);
  102.     }
  103.     for ( i = 0; i < u->dim; i++ )
  104.     if ( v_entry(u,i) < 0 || (i < u->dim-1 &&
  105.                   v_entry(u,i+1) > v_entry(u,i)) )
  106.         break;
  107.     if ( i < u->dim )
  108.     {
  109.     errmesg("svd()");
  110.     printf("# SVD sorting error\n");
  111.     }
  112.  
  113.  
  114.     /* test of long vectors */
  115.     notice("Long vectors");
  116.     x = v_resize(x,100000);
  117.     y = v_resize(y,100000);
  118.     z = v_resize(z,100000);
  119.     v_rand(x);
  120.     v_rand(y);
  121.     v_mltadd(x,y,3.0,z);
  122.     sv_mlt(1.0/3.0,z,z);
  123.     v_mltadd(z,x,-1.0/3.0,z);
  124.     v_sub(z,y,x);
  125.     if (v_norm2(x) >= MACHEPS*(x->dim)) {
  126.        errmesg("long vectors");
  127.        printf(" norm = %g\n",v_norm2(x));
  128.     }
  129.  
  130.     mem_stat_free(1);
  131.  
  132.     MEMCHK();
  133.  
  134.     /**************************************************
  135.     VEC        *x, *y, *z, *u, *v, *w;
  136.     VEC        *diag, *beta;
  137.     PERM    *pi1, *pi2, *pi3, *pivot, *blocks;
  138.     MAT        *A, *B, *C, *D, *Q, *U;
  139.     **************************************************/
  140.     V_FREE(x);        V_FREE(y);    V_FREE(z);
  141.     V_FREE(u);        V_FREE(v);    V_FREE(w);
  142.     V_FREE(diag);    V_FREE(beta);
  143.     PX_FREE(pi1);    PX_FREE(pi2);    PX_FREE(pi3);
  144.     PX_FREE(pivot);    PX_FREE(blocks);
  145.     M_FREE(A);        M_FREE(B);    M_FREE(C);
  146.     M_FREE(D);        M_FREE(Q);    M_FREE(U);
  147.     M_FREE(Atmp);    M_FREE(Qtmp);
  148.  
  149.     MEMCHK();
  150.     printf("# Finished torture test\n");
  151.     mem_info();
  152.  
  153.     return 0;
  154. }
  155.  
  156. FileDataŵtutadv†Eÿÿÿ@ò XÂ
  157. /* routines from the section 8 of tutorial.txt */
  158.  
  159. #include "matrix.h"
  160.  
  161. #define M3D_LIST    3      /* list number */
  162. #define TYPE_MAT3D  0      /* the number of a type */
  163.  
  164. /* type for 3 dimensional matrices */
  165. typedef struct {
  166.     int l,m,n;    /* actual dimensions */
  167.     int max_l, max_m, max_n;    /* maximal dimensions */
  168.     Real ***me;    /* pointer to matrix elements */
  169.                    /* we do not consider segmented memory */
  170.         Real *base, **me2d;  /* me and me2d are additional pointers 
  171.                 to base */
  172. } MAT3D;
  173.  
  174.  
  175. /* function for creating a variable of MAT3D type */
  176.  
  177. MAT3D *m3d_get(l,m,n)
  178. int l,m,n;
  179. {
  180.   MAT3D *mat;
  181.   int i,j,k;
  182.  
  183.   /* check if arguments are positive */
  184.   if (l <= 0 || m <= 0 || n <= 0)
  185.     error(E_NEG,"m3d_get");
  186.  
  187.     /* new structure */
  188.   if ((mat = NEW(MAT3D)) == (MAT3D *)NULL)
  189.     error(E_MEM,"m3d_get");
  190.   else if (mem_info_is_on()) {
  191.     /* record how many bytes is allocated */
  192.     mem_bytes_list(TYPE_MAT3D,0,sizeof(MAT3D),M3D_LIST);
  193.     /* record a new allocated variable */
  194.     mem_numvar_list(TYPE_MAT3D,1,M3D_LIST);
  195.   }
  196.  
  197.   mat->l = mat->max_l = l;
  198.   mat->m = mat->max_m = m;
  199.   mat->n = mat->max_n = n;
  200.  
  201.     /* allocate memory for 3D array */
  202.   if ((mat->base = NEW_A(l*m*n,Real)) == (Real *)NULL) 
  203.     error(E_MEM,"m3d_get");
  204.   else if (mem_info_is_on())
  205.     mem_bytes_list(TYPE_MAT3D,0,l*m*n*sizeof(Real),M3D_LIST);
  206.  
  207.     /* allocate memory for 2D pointers */
  208.   if ((mat->me2d = NEW_A(l*m,Real *)) == (Real **)NULL)
  209.     error(E_MEM,"m3d_get");
  210.   else if (mem_info_is_on())
  211.     mem_bytes_list(TYPE_MAT3D,0,l*m*sizeof(Real *),M3D_LIST);      
  212.  
  213.     /* allocate  memory for 1D pointers */
  214.   if ((mat->me = NEW_A(l,Real **)) == (Real ***)NULL)
  215.     error(E_MEM,"m3d_get");
  216.   else if (mem_info_is_on())
  217.     mem_bytes_list(TYPE_MAT3D,0,l*sizeof(Real **),M3D_LIST);
  218.  
  219.       /* pointers to 2D matrices */
  220.   for (i=0,k=0; i < l; i++)
  221.     for (j=0; j < m; j++)
  222.       mat->me2d[k++] = &mat->base[